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ABSTRACT 

For  simultaneous  non-linear  equations  the  convergence  of 
all  functional  iterative  procedures  is  dependent  upon  a  good 
initial  approximation  to  the  desired  root.  In  this  report, 
the  restriction  on  the  choice  of  the  initial  approximation 
has  been  eased  by  the  use  of  a  technique  known  as 
Parameter-Perturbation . 

Parameter-Perturbation  divides  each  problem  into  a  numbe'r  of 
subsidiary  problems,  each  of  which  is  solved  using  an  itera¬ 
tive  procedure.  Of  tfye  several  functional  iterative  procedures 
available,  Newton's  method  appears  to  be  the  most  preferable 
for  the  systems  of  equations  considered  in  this  report. 

The  solution(s)  to  the  finite-difference  equations  describing 
(a)  a  first-order  reaction  with  heat  effects  and  (b)  a 
Langmuir-Hinshe lwood  type  of  rate  mechanism,  taking  place  in  a 
spherical  catalyst  were  found  to  agree  with  the  published 
results  when  a  sufficiently  large  number  of  grid  points 
were  employed.  Due  to  the  nature  of  the  equations  describing 
the  Langmuir-Hinshe lwood  rate  mechanism,  convergence  took  more 
computation  time  than  that  used  by  the  first-order  reaction. 

Besides  the  overall  effectiveness  factor,  the  solutions  also 
give  concentration(s),  temperatures  and  local  effectiveness 


factors  as  a  function  of  radii  and  boundary  conditions. 
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I  INTRODUCTION 

Many  chemical  engineering  problems  give  rise  to  non-linear 
simultaneous  algebraic  or  transcendental  equations.  In  the 
past  a  straight  forcing  procedure  or  Newton's  method, 

Appendix  Ai ,  have  been  used  to  obtain  solutions  to  these  sets 
of  equations.  It  has  been  found  that  these  methods  have  only 
limited  applicability,  because  they  depend  upon  a  good 
initial  estimate  to  the  desired  root(s),  an  estimate  which  is 
usually  not  available.  For  real  systems  with  real  solutions, 
the  size  of  the  domain  of  convergence  is  inversely  related 
to  the  degree  and  number  of  equations. 

To  increase  the  domain  of  convergence,  Freudenstein  and 
Roth  (l)  have  suggested  an  algorithm  known  as  Parameter- 
Perturbation.  This  algorithm  essentially  consists  of 
using  a  set  of  'derived'  equations  which  are  similar  in  form 
to  the  equations  to  be  solved.  These  'derived'  equations 
have  a  number  of  arbitrarily  chosen  coefficients,  which  may 
be  'perturbed',  so  that  they  can  be  solved  by  a  standard 
iterative  procedure.  The  'derived'  equations  are  then 
linearly  deformed  into  the  'true'  equations  in  a  finite 
number  of  steps. 

The  main  purpose  of  this  report  is  to  show  how  Parameter- 
Perturbation  can  be  used  to  solve  the  large  systems  of 
simultaneous  non-linear  equations  associated  with  the 
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mathematical  model  of  a  catalyst  particle.  The  mathematical 
model  involves  dividing  the  catalyst's  radii  into  N  segments, 
with  each  segment  assumed  to  have  constant  physical  para¬ 
meters.  Non-linear  terms  arise  due  to  the  type  of  reaction 
rate  equations  employed  in  the  mathematical  model. 
Specifically  considered  in  this  report  are  a  number  of  rate 
expressions,  which  are  functions  of  both  temperature  and 
concentration.  The  technique  employed  to  solve  the 

l 

mathematical  model  involves  replacing  the  non-linear  mass 
and  energy  differential  equations  with  third-order  correct 
finite-difference  equations,  Appendix  Aii.  Solution  of 
these  finite-difference  equations  makes  use  of  the 
Parameter-Perturbation  technique . 
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II  PARAMSTNR-PHRTURBATION 

There  is  no  established  method  for  solving  sets  of  simultaneous 
non-linear  algebraic  equations.  In  many  cases,  however,  the 
type  of  solution  will  be  known  from  physical  considerat ions 
and  iterative  methods  can  be  used.  Whatever  method  is  used  to 
obtain  a  first  approximation,  successive  corrections  are  usually 
obtained  from  linear  equations.  If  the  equations  to  be  solved 
may  be  written  as: 

F  (x)  -  X  -  f  (X)  *  0 

or 

X  -  f.  (X) 

then  the  simplest  iterative  technique  is  given  by  the  recursive 
formula: 

yn-n  .  f  (x(nb 

where  the  superscript  designates  the  appropriate  iterate.  This 
particular  technique  exhibits  geometric  or  first-order  convergence. 

Another  common  iterative  technique  is  application  of  Newton’s 
method  to  function  space  and  is  given  by: 

X(n'l}  =  X(n'  -  (J-1)^  F  (X^) 

where  is  the  Jacobian  matrix  of  partial  derivatives.  This 
method  exhibits  quadratic  convergence,  and  when  it  does  converge, 
does  so  faster  than  the  previous  method. 

To  solve  a  set  of  non-linear  algebraic  equations,  assume  an 

initial  trial  vector  X^0;.  Then  using  Newton's  method  obtain 

X^n^  where  *  X^n^±  specified  tolerance,  in  n  finite 

(n) 

steps,  X  is  the  'true'  solution  vector.  However, 


' 


. 


. 


■ 


>  •  c 


('*) 


if.  the  difference  between  corresponding  elements  of  X  ^ 
and  ^  becomes  larger  as  the  number  of  iterations 

increase  then  Newton's  will  not  yield  a  solution.  In  this 
case  consider  the  algorithm,  Parameter-Perturbation. 


The  Parameter-Perturbation  procedure  is  an  algorithm  for 
determining  a  root  of  a  set  of  simultaneous  non-linear 
equations.  It  is  a  method  essentially  independent  of  the 
need  of  a  good  initial  approximation. 


To  determine  a  solution  of  a  set  of  equations  (X)  •=  0, 

i=l , 2 , 3 , . . • , n  where  X  is  an  n-dimensional  vector  with 

components  x^,  x^, . . . ,x^  and  (X)  is  of  the  form 
m . 

F1<2)  =  ^  piAk<^ 

k=o 


where  the  P.,  are 
lk 


of  equations  the 


parameters,  first  consider  another  set 
'derived'  equations  G^°^(X)  =  0  ,  i=l,2,.. 


These  'derived'  equations  may  be  any  set,  with  one  known 

root,  which  belongs  to  the  same  family  as  F^ (x)  ,  that  is 

m . 


g(°)(x)  = 

where 


P .  ^  ,  except 


for  i=k. 


(  o  \ 

The  'derived'  equations,  G^  '  ^ (X)  =0,  are  deformed  into 
the  equations  F-^  (X)  =0  by  means  of  a  finite  number,  N,  of 


b 
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successive  small  increments  in  the  parameters.  Define  N  sets 

of  equations 

m . 

1 

Gi(j)(x)  Qik(j)?ik(X)  .  j  =  i,2,...,N 

k  =  o 

such  that 


Gi(N)(X)  =  F.(X)  , 


Q 


ik 


( J  )_ 


=  Q  -i 

xik 


(o) 


+  (P..  -Q., 
v  ik  ^ik 


(o) 


)i 

N 


i 


The  method  does  not  require  the  steps  to  he  of  equal  size. 
The  desired  root  is  obtained  by  solving  the  N  sets  of 
equations  using  Newton's  method.  The  algorithm  converges  to 
a  root  of  F.(X)  =0  if  the  functions  F.(X)  and  G.^°^(X)  are 
such  that : 

(1)  G^(t,x)=0  is  continuous  for  0-  t  ~N 

(2)  Xf(-t)  is  continuous  for  0-  t  -N  . 


To  decrease  the  computation  time  it  is  important  to  choose 
the  maximum  increment  for  which  Newton's  method  will 
converge.  If  upon  choice  of  an  increment  size  no 
convergence  occurs,  the  increment  is  halved  for  the  next 
trial.  If  convergence  occurs  for  two  successive  trials, 
then  for  the  next  trial  increment  is  doubled.  An 
example  of  this  procedure  is  given  in  Table  1,  page  7. 
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An  example  problem  is  worked  out  in  Appendix  B,  using 
Parameter-Perturbation  with  Newton's  method.  The  results 
show  that : 

(1)  Newton's  method  converges  faster  than  the  function 
iteration  procedure. 

(2)  The  solution  is  obtained  using  any  physically  feasible 
trial  vector. 
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TABLE  1 

Convergence 

Proc 

Trial 

Number 

'Total ' 
Increment 

Incr 

S 

1 

1.000 

1. 

2 

.500 

• 

3 

.250 

• 

4 

.125 

0 

5 

.250 

0 

6 

.500 

0 

7 

1.000 

0 

8 

-750 

0 

9 

1.000 

0 

10 

.875 

• 

11 

1.000 

• 

t  Comment 


No  convergence,  reduce 
increment  by  one-half 

No  convergence,  reduce 
increment  by  one-half 

No  convergence,  reduce 
increment  by  one-half 

t 

Convergence,  retain 
same  increment 

Convergence,  double 
the  increment 

Convergence,  double 
the  increment 

No  convergence,  reduce 
increment  by  one-half 

Convergence,  retain 
same  increment 

No  convergence,  reduce 
increment  by  one-half 

Convergence,  retain 
same  increment 

Convergence,  solution 
obtained 


edur 

emen 

i  ze 

000 

500 

250 

125 

125 

250 

500 

250 

250 

125 

125 
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III  MATHEMATICAL  MODEL  OF  A  CATALYST  PARTICLE 

Consider  a  catalyst  pellet  with  radial  symmetry,  and  let  us 
focus  our  attention  on  a  shell  of  thickness  dr  and  radius  r. 
Assume  that  the  complicated  diffusion  phenomena  within  the 
porous  structure  obeys  Pick's  law  which  assumes  that  the 
diffusivity  coefficient  is  constant  for  a  given  pressure, 
temperature,  and  pair  of  substances,  and  essentially 
independent  of  concentration. 


Reactants  are  transported  to  and  from  the  shell  by  diffusion 
and  consumed  within  the  differential  shell  by  reaction.  At 
steady-state,  a  mass  balance  for  some  component  A  is; 

(5-1) 

(Rate  of  diffusion  of  A  inward  at  r=r+dr)  -  (Rate  of 
diffusion  of  A  inward  at  r=r)  +  (Generation  of  A  due  to 
reaction  in  the  shell)  =  0 


In  differential  form  this  becomes; 


(3-2) 


4il(r+dr)2D  (dC.  +d2C .  dr  )-4Tlr2D  dC.-4nr2drS  f  .  (C  .  ,  C„ ,  .  .  ,  T  )=0 
v  /  Y\  a  A  '  r  A  v  lv  A’  B’  ’  ' 


dr  dr‘ 


dr 


2 

where  4Tir  is  the  inner  superficial  area  of  the  spherical 

shell,  D  is  the  diffusivity,  dC.  is  the  concentration 
’  r  J  ’  _ A 

dr 

gradient  at  the  shell  surface,  and  Sv  is  the  pore  surface, 

2  t 

(cnf/cnr  of  the  volume  of  the  porous  structure). 


■ 

— 


(9) 


The  above  equation  (3-2)  reduces  to; 


(3-3) 


A 


4-  2  dC 


A 


dr 


2 


r  dr 


I) 


r 


(3-'i ) 


A’  B’  *  *  *  ’ 


A  mass  balance  on  component  (i),  thus  yields; 


(3-5) 


+  2  dC . 

l 


r  dr 


Consider  the  same  catalyst  pellet  as  that  used  in  deriving 
equation  (3-5)  and  assume  heat  transfer  by  conduction 
within  the  shell  obeys  Fourier's  law,  which  assumes  that 
thermal  conductivity  k  is  a  function  of  the  molecular 
state  of  the  medium. 

Heat  is  either  consumed  or  produced  within  the  shell  due 
to  the  chemical  reaction(s),  and  transported  to  or  from  the 
shell  by  conduction.  At  steady-state  an  energy  balance  on 
the  shell  becomes; 


(3-6) 


(Rate  of  conduction  inward  at  r=r+dr)  -  (Rate  of  conduction 
inward  at  r=r)  -  (Rate  of  heat  generation  in  shell)  =  0 
In  differential  form  equation  (3-6)  becomes; 


(3-7) 


M(r  +  dr)2K(dT  +  .  d2T  dr)  -  4Tlr2K  dT  -  4~ir2dr  f  ( CA ,  CB 


T)=0 


dr 


dr 


l 


-  m 


/ 


' 


. 


. 


. 


/ 


(10) 


where  K  is  the  effective  thermal  conductivity  coefficient, 

dT  is  the  thermal  gradient  at  the  shell  surface, 
dr 

The  above  equation  (3-7)  reduces  to; 


(3-8) 


d2T  2  dT  \  Hifi(CA,CB,...,T) 

2  ~  " — 

dr  r  dr  l  K 

Reasoning  similar  to  that  used,  in  obtaining  equations  (3-5) 
and  (3-8)  for  a  spherical  shell,  enables  one  to  derive  the 
corresponding  equations  for  cylindrical  and  slab  geometry. 


For  cylindrical  geometry  the  equations  are; 


2 

d  C. 

l 

+ 

1  dC. 

l 

fi (CA’CB’ * ' * 

,T) 

,  2 
dr 

r  dr 

D 

r 

and 

2 

d  T 

+ 

1  dT 

\~  H .  f  .  (C.  ,0^, 

.  .  .  ,T) 

%  2 

dr 

r  dr 

i  K 

where 

i 

refers 

to 

components  A, 

6 ,  ... 

For  a 

one-dimensional  slab  the 

equations 

d2C  . 

l 

fi<CA’ 

CB- 

• • .,T) 

CM 

XI 

"d 

D 

X 

(3-9) 


(3-10) 


(3-H) 


and 


’ 


■ 


/ 


(11) 


(3-12) 


9 

d  T 


Hifi(CA’CB’ ' • ‘ ’T) 


dx  i  K 

These  equations  have  the  following  boundary  conditions; 


At  r  =  R 


C.  =  C. 

1-  is 

T  =  T 


At  r  =  0 


dC .  =  0 

i 

dr 

dT.  =  0 
dr 


At  x  =  L 


C.  =  C. 
r  is 

T  =  T 


At  x  =  0 


dC.  =  0 

_i 

dx 

dT  =  0 
dx 


(3-13) 


The  grid  system  employed  for  a  spherical  catalyst  is  of  the 
following  form; 


\ 

i 

(l+2n)R 

2N 

i 

i 

/ 

/ 

/ 

/ 

/ 


"  -  ■'  ja  -m 


. 


. 


(12) 


where  N  is  the  total  number  of  grid  points 

n  position  of  grid  points,  calculated  from  the  centre. 

The  finite-difference  approximation  for  the  space  derivatives 
are  : 


(3-15) 


9 

d“T 


dr 


2 


T  ,  -  2  T  +  T  , 
_  n+1  n  n-1 

2 


n  =  1,2,... , N-1 


(3-16) 


dT  2T,+3T-6T,+T0  ,  n 

=  n+1 _  n _ n-1 _ n-2  n  =  1,2,...,N-1 


dr 


6£r 


for  n=N,  outermost  grid  point,  the  space  derivative  are; 


2 

d  T 


dr‘ 


— TN— 2  +  10  TN_!  -  25  TN  +  16  TS 

5Ar2 


(3-17) 


dT 

dr 


3  Tn_2  -  20  TN-1  -  15  Tn  +  32  Tg 

30ir 


(3-18) 


Similarity  for  the  mass  balance  we  have; 


A. 

17' = 


!.  .  -  2  C.  +  C. 

l , n+1  i,n  l , n-1 

A  r2 


(3-19) 


. 


EP 


> 

lb 

■ 

■ 

. 


. 


* 

*■ 
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dC.  2  C.  .  +  3  C.  -  6  C.  ,+C.  . 

i  =  i ,  n+1  l ,  n _ i,n-l _ l ,  n-2 


dr 


6Ar 


2 

d^C  . 

l 


-  C.  M  n  +  10  c.  XT  „  -  25  c.  r  +  16  c.  c 

i,N-2  i,N-l  i,N  i,S 


dr' 


5f.\r‘ 


dC.  3  C.  ,T  0  -  20  C.  ,T  ,  -  15  C.  ,T  +  32  C.  c 

l  =  l  , N-2  _ l  ,  N - 1 _  l  ,N  l  ,  S 


dr 


30Ar 


(3-20) 


(3-21) 


(3-22) 


where  n  =  1,2,...,N-1 

i  =  components  of  the  system 


Once  the  set  of  finite-difference  equations,  approximating 
the  mathematical  model's  differential  mass  and  energy 
balances,  have  been  solved  for  the  concentration  and 
temperatures,  it  is  possible  to  obtain  the  local  and 
overall  effectiveness  factors  for  the  catalyst  particle. 


The  local  effectiveness  factor 


is  defined  as ; 


Rate  of  reaction  at  r  =  r 


Rate  of  reaction  at  r  =  R 


or 


(3-23) 


The  overall  effectiveness  factor 


9 


is  defined  as ; 


' 


. 


i  ;  l©  rn  .toque  * 


' 

' 


' 


.  .  ,  T  )  dr 

T  ) 

’  S ; 


(3-24) 


It  was  decided  to  integrate  equation  (3-24)  piecewise, 
fitting  exactly  through  the  points  a  fourth  degree 
polynomial  to  each  segment.  The  derivation  is  given  in 
Appendix  Aiii. 


The  rate  constant(s)  are  considered  to  he  a  function  of  the 
absolute  temperature  and  the  activation  energy  associated 
with  the  reaction,  i.e. 

(3-25) 

(  5-  -  —  ) 

k  =  k  e  RT  RT 
v  vs  g  s  g 


The  derivation  of  equation  (3-25)  is  given  in  Appendix  Aiv. 
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IV  DETAILS  OF  SOLUTION 

(a)  First  Order  Irreversible  Reaction  With  Heat  Effects 

To  show  the  applicability  of  Parameter-Perturbation  in 
finding  the  solution  to  systems  of  simultaneous  non-linear 
equations  we  have  chosen  as  our  first  system  a  mathematical 
model  describing  an  irreversible,  first  order  reaction  with 
heat  generation.  The  equations  describing  this  first  order 
reaction  are  as  follows:  , 

Stoichiometry  (4-l) 

X  - Products 


Rate 


r  =  k  e 
X  vs 


E 


E 


R  T 
g  s 


) 

E  T.  . 


(h-2) 


The  solutions  to  the  mathematical  model  of  an  irreversible 
first  order  reaction  with  heat  generation  have  been 
presented  by  a  number  of  people,  some  of  whom  are;  Weisz  and 
Hicks  (3),  Tinkler  and  Metzner  (4),  and  Carberry  (5).  However 
none  of  the  solutions  presented  up  to  now  have  made  any  use 
of  the  finite-difference  technique  employed  in  this  report. 

In  order  to  present  the  solution  in  a  standard  form  it  is 
necessary  to  define  the  following  dimensionless  parameters: 
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Employing 
with  the 
to  (4-5) 
was  used 
catalyst 
parameter 

The  solut 
irreversi 
versus 
represent 
and  -.20, 
by  Weisz . 
Table  2, 

Using  lira 
70,000,  i 
could  be 


('*-3) 


(**-'*) 


S 


H  D 


xs  r 


T  K 
s 
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the  mathematical  model  presented  in  Section  III 
specific  modifications  as  given  by  equations  (4-1) 
inclusive  the  Parameter-Perturbation  algorithm 
to  obtain  temperatures  and  concentrations  in  the 
particle  as  a  function  of  the  dimensionless 
s  ^ and  \j  . 

ions  to  the  mathematical  model  describing  an 
ble  first  order  reaction  are  presented  as 
,  (j>  and  j-)  .  Figure  1,  page  18  shows  the  graphical 

ation  of  these  results  for  f\  =  30 ,  R  =  0.0,  .05 

also  shown  are  the  corresponding  values  obtained 
and  Hicks.  Numerical  results  are  compared  in 
page  34 . 

ited  numerical  accuracy,  core  storage  limited  to 
t  was  observed  that  no  numerical  convergence 
obtained  for  cases  where  A  '  +.15.  No  problems 
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of  convergence  were  encountered  with  endothermic  first 
order  reactions.  Refer  to  Table  3,  page  37. 

With  the  transformation  used  by  Weisz  and  Hicks  only  the 

overall  effectiveness  factor  can  be  obtained, but  with 

the  finite-difference  approach  concentration , temperature 

and  effectiveness  factor  as  a  function  of  the  reduced 

radii  can  be  solved  for. Example  plots  of  reduced  concentration 

effectiveness  factor  and  temperature  versus  reduced  radii 

t 

are  given  in  Figures  2,  3  and  4  respectively.  A  more 
complete  tabulation  of  results  is  available  from  the 
Chemical  Engineering  Department.  Figures  2  ,  3  and  4  show 
that  the  major  change  in  concentration  and  temperature 
takes  place  in  the  outer  one-third  of  the  spherical  catalyst. 
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FIGURE  1 

Comparison  of  Effectiveness  Factors  in  Terms  of  Observable 

Parameters 
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Reduced  Concentration  versus  Reduced  Radii, First  Order 
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Local  Effectiveness  Factor  versus  Reduced  Radii, First  Order 


FIGURE  3 
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Reduced  Temperature  versus  Reduced  Radii,  First  Order 
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(b)  Lan ymuir-iiinshelwood  reaction  With  Heat  Effects 


for  our  second  example  of  the  applicability  of 
Parameter-Perturbation  we  propose  to  find  the  solution(s) 
to  the  mathematical  model  of  a  reaction  taking  place  in 
a  spherical  catalyst,  where  the  rate  expression  is  of 
the  Langmuir-Hinshelwood  type. 


The  general  chemical  equation  describing  reactions 
under  consideration  are; 

A  +■  bi->  -+-...  •+•  - ►  xX  ■+  yY  +  ... 


(4-6) 


The  rate  equation  is  taken  to  be, 

r  =  kpA/UtXAPA+  Zxipi) 

i 

where  i  denotes  any  reaction  product  or  reactant  other  ' 


(4-7) 


than  A. 


Effectiveness  factors  for  an  isothermal  catalyst  have 
recently  been  presented  by  Roberts  and  Satterfield  for 
a  catalyst  of  slab  geometry  (6).  Their  method  of 
solution  can  be  summarized  as  follows.  First  assume:  a  slab 
geometry;  ideal  gas  behaviour;  an  effective  diffusivity  for 
each  species  which  is  constant  but  not  necessarily  equal  to 


each  other. 
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A  material  balance  on  component  A,  over  a  differential 


thickness 

within 

the  catalyst, 

DAd2C 

A  = 

da(  1 

)  d2pA  =  r 

CM 

X 

RT 

dx2 

A  similar 

balance 

2 

on  any  other 

Di(  —  ) 

d  P. 

-V.  r 

l  = 

xi 

RT 


dx‘ 


(4-8) 


(4-9) 


where  is  the  stoichiometric  coefficient. 


Combine  equations  (4-8)  and  (4-9)  and  integrate  subject 
to  the  boundary  conditions  p  =p.  ;  p.=p.  at  x=0 

A  A, S  1  1  ,  S 


and 

dFA 

=  =  o 

at  x  =  Li 

dx 

dx 

(4-10) 


to  yield  an  expression  for  p.  in  terms  of  p. . 

1  A 

Substituting  the  resulting  expression  into  equation 
(4-7)  gives 

(4-11) 


r  =  kp 


(1  +  Pa(Ka-Da  L  (V^/Di))+  t  Ki(Pi,st(fA,sW1)i>» 


Define  w  =  1  +  %  K±  (PjL  (  s  +  ( PA ;  s  '£VDi  ) ) 


(4-12) 


Note:  Only  positive  values  of  w  can  be  handled. 
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Define  k'  =  k/w 


(2k) 

(k-13) 


and 


K  = 
o 


(vDA^KiY/Di»/’' 


(4-14) 


Using  these  definitions,  equation  (4-ll)  reduces  to 
r  =  k’pA/(i+KoPA) 


(k-15) 


Equation  (4-15)  way  he  substistuted  into  equation  (4-8) 
and  integrated  to  give 


d(KoPA) 

d(x/L) 


"  'I2  '(*m^Ko*pA_PA  0)  ~  ln(  - ~  ))“ 

II  °  A  A,0  i+K  p 

A  > 


(4-16) 


where  0^  =  L 


/  k  1 ET  n.  2 
V  Da  } 


In  terms  of  present  nomenclature  the  effectiveness 
factor  is  then  defined  as 


\ 


C?  i+K  p, 
^2  ^  o*A , 

0  '  K  p. 

o^A,s 


1+K  p.  i 

)(K  (p.  -p.  )-ln( - 2-Ai®.  ))2 

°  A,s  "A, o  ^  ' 

o^A ,  o 


( k-17 ) 


where  p.  is  determined  by  numerically  integrating 

j  O 

equation  (4-l6) 


where 

p  o  is  the  partial  pressure  at  the  center  of  the 
catalyst  particle 
Kq  overall  adsorption  constant 
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k  rate  constant 

k'  modified  rate  constant 

PA  g  partial  pressure  at  the  surface 

Employing  the  mathematical  model  as  defined  in  Section  III, 
which  differs  from  that  used  hy  Roberts  and  Satterfield 
because  the  former  is  for  a  spherical  catalyst  pellet  with 
heat  generation.  For  spherical  coordinates  the  solution 
is  as  follows; 


Mass  balance 

(4-18) 

Dr  ,  d2PA  2  dPA,  k'PA 

R  T  dr  r  dr  1+K  p. 

g  o5_A 


Similarily  the  energy  balance  is 


(4-19) 


K  ( 


d2T 


dr 


2  dT  AH.klp 
2  +  ~  =  -  1  -  - 


r  dr 


1+K  p. 
o^A 


where 


(  -  • 
/kv  '  R  T 

(*)  e  g  s 


E 

R-  T 


) 


Solutions  using  Parameter-perturbation  were  obtained  for 
the  system  C,  CO^  previously  investigated  by  Austin  and  Walker 
(7).  However  Kq  was  held  constant  for  the  solutions  obtained, 
the  error  introduced  by  this  simplification  is  not  known. 


Results  of  Cj  versus  (j)  ,  R  and  p^gK^ ,  obtained  by  using 
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Parameter-Perturbation  to  solve  the  sets  of  finite- 
difference  equations  representing  both  spherical  and 
slab  geometry  are  presented  for  the  system  carbon, 
carbon-dioxide  in  Figure  5  page  28  . 

The  approach  employed  by  Roberts  and  Satterfield  applies 
only  to  slab  geometry,  and  for  reactions  without  any  heat 
generation.  As  noted  from  equation  (4-9)  only  the  overall 
effectiveness  factor  can  be  evaluated.  However  with  the 

i 

finite-difference  approach  any  type  of  geometry  can  be 
used,  as  well  as  reactions  with  heat  generation. 
Concentration,  effectiveness  factor  and  temperature  are 
solved  for  as  functions  of  the  reduced  radii. 

Example  plots  of  reduced  concentration,  effectiveness  factor 
and  temperature  versus  reduced  radii  are  given  in  Figures 
6,  7  and  8  pages  29  ,30  and  31  respectively,  for  p.  K  =  -.964 
A  more  detailed  tabulation  of  results  are  available  from  the 
Chemical  Engineering  Department. 


From  the  calculations  carried  out  using  both  spherical 

and  slab  geometry  for  pA  K  =  -.964  it  was  observed  that: 

J  ^As  o 

(a)  Effect  of  heat  generation  on  overall  effectiveness 


o, 


factor  is  neglible  for  a  surface  temperature  of  1273  K, 

assuming  K  to  be  constant, 
o 


I 


' 

j.i 

• „  * 
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■ 

■  . 
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(b)  The  effectiveness  factors  calculated  for  p.  K  =  -.96*1 
are  much  lower  than  the  corresponding  effectiveness 
factors  for  a  first  order  reaction  with  identical 
physical  parameters. 

(c)  Plots  of  overall  effectiveness  factors  as  a  function 

of  the  parameters  j'  and  p.  K  for  the  cases  solved 
^  x’  y  rAs  o 

hy  both  Roberts  and  myself  coincide,  but  if  we  make 

2  g 

use  of  the  transformation  L  =  R“/9  proposed  by 
R.  Aris  (9)  to  convert  from  linear  to  spherical 
coordinates  the  results  obtained  from  linear 
coordinates  do  not  coincide  with  my  results  for 
the  spherical  system.  Thus  one  is  led  to  the 
conclusion  that  for  this  particular  rate 
mechanism  and  pA  K  =  -.964  the  transformation  does 
not  hold. 

From  Figures  6,  7  and  8  pages  29,  30  and  31  respectively 
the  following  is  observed: 

(a)  Concentration,  effectiveness  factor  and  temperature 
gradients  are  steepest  near  the  surface. 

(b)  Concentration  near  the  centre  of  the  particle  is 
high  in  comparison  to  the  corresponding  effectiveness 


factor 


* 


. 


FIGURE  5 


(28) 


Comparison  of  Effectiveness  Factors  in  Terms  of  Observable 

Parameters 
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Partial  Pressure  versus  Reduced  Radii  and  Length , Langmuir-Hinshe lwood 


FIGURE  6 
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and 
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Local  Effectiveness  Factor:  versus  Reduced  Radii  and  Length,  Langmuir- 
Hinshelwood 


FIGURE  7 


(•30) 


- 


' 


• 

Reduced  Temperature  versus  Reduced  Radii  and  Length,  Langmuir-Hinshe lwood 


FIGURE  8 


(•51) 


(32) 


V  RESULTS 
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of  the  catalyst  consist  of  the  known  surface 
concentration ( s )  and  temperature. 

(2)  Thermal  conductivity  and  diffusion  coefficients 
can  be  fed  into  the  program  as  a  function  of  the 
radii,  if  this  function  is  known. 

(3)  Concentration,  temperature  and  effectiveness  factor 
are  calculated  as  functions  of  the  reduced  radii. 

To  minimize  round-off  error  the  finite-difference 
equations  are  multiplied  by  arbitrary  constants,  such 
that  all  the  equations  have  parameters  which  are  of  the 
same  order  of  magnitude.  For  the  systems  considered  in 
this  report  we  chose  1  as  the  order  of  magnitude  with 
which  to  work.  Convergence  to  an  acceptable  solution  for 
the  individual  grid  system  being  considered  is  assumed 
whenever  the  residual  of  each  equation  is  less  than  . 0005. 
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The  accuracy  with  which  the  system  of  finite-difference 
equations  are  solved  by  the  computer  is  dependent  only 
upon  the  accuracy  with  which  the  ’true’  set  of  equations 
are  solved,  and  not  upon  the  number  of  'deformed'  sets  of 
equations  used  to  arrive  at  the  'true'  set  of  equations. 


Table  2  is  a  s 
by  the  Paramet 
are  compared  t 
for  =-.20 
solutions  for 
are  all  for  a 
numerical  conv 
of  center  line 
grid  points,  f 
Hinshelwood  ra 


ummary  of  some  of  the  solutions  obtained 
er-Perturbation  algorithm.  These  results 
o  solutions  obtained  by  Weisz  and  Hicks 
and  jO  =  +.05  and  to  the  analytical 

J 

the  case  where  R  =  0,  these  comparisons 

/ 

first  order  rate  mechanism.  To  show  that 
ergence  has  been  obtained,  a  comparison 
concentrations  obtained  using  19  and  23 
or  both  First  Order  and  Langmuir- 
te  mechanisms  are  shown. 
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TABLE 


2 


Error  Calculations  for  Effectiveness  Factor 


,(3'0 


Rate  Mechanism 

$ 

n<1) 

(2) 

q 

^Error 

C3) 

First  Order 

1.778 

0.00 

.889 

.888 

.11 

(Spherical  - 

6.537 

0.00 

.654 

.652 

.31 

Coordinates ) 

12.356 

0.00 

.475 

.473 

.42 

15.587 

0.00 

.410 

.408 

-50 

2.108 

.05 

1.054 

1.054 

.00 

• 

6.005 

.05 

1.001 

1.000 

.  10 

1.275 

-.20 

.637 

.635 

.31 

6.746 

-.20 

.259 

.257 

.78 

i 

Rate  Mechanism 

D 

R 

k 

c 

C  /5) 

Numerical 

X 

V 

vs 

0 

Error  (6) 

First  Order 

.005 

0.00 

.010 

.731 

.731 

0.0 

(Spherical  - 

.005 

0.00 

.  190 

.026 

.026 

0.0 

Coordinates) 

.005 

-.20 

.010 

.849 

.849 

0.0 

.005 

-.20 

.  190 

.552 

.552 

0.0 

.005 

+ 

• 

0 

.010 

.640 

.640 

0.0 

Rate  Mechanism 

D 

X 

k 

vs 

5  X  10~7 

Pa  K 

As  0 

P  (4) 

a  4 

P  (5) 
A  4 

Numerical 

Error 

Langmuir- 

.025 

-.964 

.354 

.354 

0.0 

Hinshe lwood 

(Spherical  - 
Coordinates ) 

.225 

5  X  10-7 

-.964 

.602 

.602 

0.0 

where : 

(1)  Overall  effectiveness  factor  calculated  using 
Parameter-Perturbation  and  23  grid  points. 

(2)  Overall  effectiveness  factor  as  calculated  by 
Weisz  and  Hicks  (  P  =  -.20  and  +.05)  or  by 
analytical  means  (  \~]  =  0.0). 
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(35) 


(3) 

(4) 

(5) 

(6) 


%  Error 


Calc  u late  d 


A  -  H 


Published 


Published  C 


Centerline  concentration  calculated  using 
23  grid  points. 

Centerline  concentration  calculated  using 
19  grid  points. 

Numerical  error  =  difference  between  results 

obtained  by  using  23  and  19 
grid  points. 


X  100 


From  the  brief  error  analysis,  Table  2,  it  is  obvious  that 
the  number  of  grid  points  required  for  a  finite-difference 
representation  of  highly  non-linear  differential  equations 
is  high.  Another  source  of  error,  aside  from  the  finite- 
difference  approximations,  is  the  numerical  integration 
procedure  required  to  obtain  an  overall  effectiveness 
factor.  However,  since  integration  is  a  smoothing  process, 
errors  generated  are  usually  small. 

Newton’s  method  of  convergence  was  found  to  be  adequate 
by  itself  whenever  only  endothermic  heat  effects  were 
considered.  Exothermic  systems  however  presented  a  much 
more  difficult  problem.  Using  limited  numerical  accuracy 
it  was  not  possible  to  solve  any  problems  for  which 

h  *  !’]  <  .15  . 

/ 


. 


' 


Without  exception  numerical  solutions  to  problems  where 

A  0  > 

U  I  /. 15  ,  the  residuals  of  the  finite-difference 
equations  could  not  he  made  to  approach  zero,  thus 
making  it  impossible  to  obtain  a  solution,  reference 
Table  3,  page  37  .  The  reason  for  this  appears  to  be  the 
limited  numerical  accuracy  of  the  computer,  total  core 
storage  of  70,000,  and  not  a  fault  of  the  method  used. 

An  inspection  of  the  finite-difference  matrix  revealed 
that  the  non-linear  term  is  of  an  order  of  one-hundred  , 
of  the  magnitude  of  the  linear  terms,  thus  to  handle 
this  problem  requires  a  large  computer  pocessing 
almost  unlimited  numerical  accuracy.  This  appears 
to  be  the  major  limit  of  the  finite-difference  approach. 

The  application  of  the  Parameter-Perturbation  technique 
to  orders  other  than  those  given  in  this  report  follows 
the  same  general  procedure  as  those  solved  here, 
reference  Appendix  C,  this  is  an  advantage  of  this 
method  of  solution  over  methods  presented  in  the 
literature.  Furthermore  the  method  outlined  in  Appendix  C 
applies  to  systems  which  have  not  been  solved  by  any 
other  approach. 

Thus  the  Parameter-Perturbation  technique  with  Newton’s 
method  is  a  general  method  of  solution  to  the  various 
proposed  rate  expressions.  The  calculation  time  required 
is  of  the  order  of  20  seconds,  I.B.M.  7020,  to  obtain 
a  complete  solution  to  any  particular  rate  expression. 


. 


. 


.  t 


i 


. 


. 

. 


■ 


(■37 ) 


TABLE 

3  Summary 

Of  Computer 

Output 

For 

• 

00 

ci 

CN 

T 

N 

R1 

Rn  r 

1 

rn 

No. 

of  Iteration 

— 

C 

c 

T 

T 

s 

s 

s 

s 

1.757 

1.016 

.824 

.995 

1.066 

.001 

.  686 

.000 

1 

1.03*i 

.996 

1,002 

1.000 

.323 

.004 

.054 

.000 

2 

.616 

.984 

1.034 

1.001 

.015 

.002 

.005 

.  000 

3 

.594 

.983 

1.033 

1.001 

.000 

'  .002 

.076 

.000 

4 

.644 

.984 

1.023 

1.001 

.006 

.002 

.343 

.000 

5 

.643 

.984 

1.023 

1.001 

.001 

.002 

.325 

.000 

/ 

6 

.632 

.984 

1.025 

1.001 

.006 

.002 

.229 

.000 

7 

.632 

.984 

1.025 

1.001 

.006 

.002 

.232 

.000 

8 

.632 

.984 

1.025 

1.001 

.006 

.002 

.  241 

.000 

9 

EAC  = 

24000. 

DKV  = 

.01 

DEX  = 

.005 

BALP  = 

.0005 

1! 

% 

10 

DELH  =  7980 


' 


' 

to  i  eso.t 


• 

. 


VI 


CONCLUSIONS 
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Finite-difference  equations  representing  a  mathematical 
model  yield  complete  solutions  for  each  grid  point,  thus 
giving  a  much  more  detailed  description  of  the  proposed 
model.  This  feature  is  of  utmost  importance  when  design 
of  equipment  is  based  upon  solutions  to  the  mathematical 
model.  Furthermore  the  use  of • Parameter-Perturbation 
enables  one  to  solve  large  sets  of  finite -difference 
equations  within  a  reasonable  length  of  time. 

In  more  detail  the  use  of  finite-difference  equations  and 
Parameter-Perturbation  to  solve  these  equations  representing 
a  boundary  value  problem  is  to  be  recommended  because: 

(1)  The  algorithm  for  obtaining  the  finite-difference 
equations  approximating  a  boundary  value  problem  and 
solution  of  these  equations  using  Parameter-Perturbation 
is  a  general  one. 

(2)  Transformation  of  variables  or  use  of  dimensionless 
parameters  is  not  necessary. 

(3)  Only  measurable  physical  quantities  are  required  to 
define  the  problem. 

■X* 

Note  Numerical  accuracy  of  the  calculations  carried  out  by 
the  computer  limits  the  type  of  problems  which  can  be 


solved . 


■  • 


' 


. 


-  * 

. 
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VII  NOMENCLATURE 


A 

reactant,  Langmuir-Hinshe 1 wood  rate 
expr e  s  s i on 

A. 

matrix  expressing  the  quantities 
rm  m  =  1 ’ 2 ’ *  *  *  » 5  ’  p  =  0 , 2 , 3 , 4 

B 

reactant,  Langmuir-Hinshe lwood  rate 
expression 

C 

concentration  of  reactant  species 

C 

X 

concentration  of  reactant  X 
(gm  moles/cc) 

D 

diffusion  coefficient  (cm  /sec) 

E 

activation  energy  (calories) 

F . (x  .  ) 
i  J 

function  of  variable  x;  i  =  l,2,...,n, 
j  =  1,2,  ...  ,n 

F  f  .  ( x  .  )  ,  F" .  (x . ) ,  .  . . 
i v  J  7  i v  J 

derivatives  of  the  function  F.(x.) 

G  . 

l 

incremental  value  used  in  parameter- 
perturbation 

K 

thermal  conductivity  of  porous 
catalyst  (cal. /sec  cm  K) 

ka 

adsorption  constant  for  reactant 

A  (atm  * ) 

K. 

1 

adsorption  constant  for  i^*1  species 
in  Langmuir-Hinshe lwood  rate 
expression  ^^-L 

K 

0 

defined  as  (Ka  -  Da  L  (K_.Vi/D.)/W 

i 

L 

length  dimension,  cartesian  coordinates 

L(l),  L(X),  ... 

linear  operator,  used  in  derivation  of 
integration  formula 

P 


partial  pressure,  (atm) 


■ 


... 


• 

'  — 


■ 


. 


, 


(40) 


jk 


^jk 


0n 

*  Jk 


R 


R 


g 


T 

W 


w. 

1 


b 

f(x) 

f  (c) 

h  . 

J 

k 


k  ’ 

k  (T) 
s  ' 


kv(T) 


matrix  expressing  the  'true'  equations 
parameters 

matrix  expressing  the  'derived' 
equations  parameters 

matrix  expressing  the  'perturbed' 
equations  parameters 

radius  of  spherical  catalyst  or  radii 
of  cylindr ial  catalyst 

universal  gas  constant 

2  3 

pore  surface  cm  /cnr  of  volume  of  the 
porous  structure 

i 

temperature  in  degrees  -  Kelvin 

maximum  incremental  value,  equal  to  one 

constant  coefficients  of  the  polynomial 
integration  formula 

product  in  Langmuir-Hinshe lwood  type  of 
rate  expression,  otherwise  the  reactant 

product  in  Langmuir-Hinshe lwood  type  of 
rate  expression,  otherwise  the  reactant 

p 

vector,  made  up  of  R  ,  p  =  1,2, 3,4, 5 

P 

any  general  function 

rate  mechanism  as  a  function  of  cone. 

residual  of  (x.  -  r.) 

x  J  J 

reaction  rate  constant,  (g .moles/cc . sec ) 
modified  rate  constant,  (g. moles  atm/cc.sec) 

rate  constant  based  on  surface  area  of 

catalyst  ,  ,  //  2  x 

J  (gin.  moles/ (cm  sec; 

rate  constant  based  on  volume  of  the 
catalys  t 


■ 


' 


('ll) 


m  . 


r 


r  . 
J 


X 

X  j  ( 1 )  ’  X  j  (  2  )  >  •  •  • 

H 

dc  or  dkD 
dr  dx 

,2  A2 

d  c  or  d  c 

T~2  ,  2 

dr  dx 


order  of  reaction  0-^.m^.l 

radii  at  position  r 
t  h  t  Ji 

j  root  of  the  i  equation,  Newton's 
method 

dimension  in  cartesian  coordinates 

first,  second,  etc.  approximation  to 

r.,  Newton's  method 
J 

heat  of  reaction  (cal/g  mole°K) 
concentration  gradient 


gradient  of  concentration  gradient 


dT  or  dT 
dr  dx 


temperature  gradient 


2  2 
d  T  or  d  T 

,2  ,2 

dr  dx 


gradient  of  temperature  gradient 


h 

0 


defined  as  C.  AH.D  E 

is  l  r 


T  2K  R 
s  g 


defined  as  C  AH.D 

s  l  r 

T  K 
s 


defined  as 


defined  as 


E 


R  T 
g  s 


/  k  C  . 

/  vs  l s 

D 


J 


* 

- 

■ 


. 


- 


■ 


■ 


* 


('i2) 


i 

j 


defined  as  (R  /D  C.  )  .  (observed 

r  is 

reaction  rate/gross  catalyst  volume) 

defined  as  the  effectiveness  factor 

defined  as  either  reaction  rate  at 
position  n  or  stoichiometric 
coefficients 

index  denoting  species 

index  denoting  position  of  an  element 
of  a  matrix 


k 


index  denoting  position  of  an  element 
of  a  matrix  ' 


n 


denotes  position,  or  number  of  the 
iteration 


s 


denotes  surface 


. 


■ 


(^3 ) 
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IX  APPENDIX 
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A 


DERIVATIONS 


(i)  Newton’s  Method 

Suppose  the  following  equations  are  to  be  solved 


(— 1 

X 

1 —> 

x2  ’ 

•  •  •  J 

y. 

p 

11 

0 

F2(x1, 

x2  ’ 

•  *  •  j 

X  )  = 

n ' 

0 

Fn(xl’ 

X 

fO 

.  .  .  , 

x  )  = 
n ' 

0 

Let  r  . 
J 

,  3  = 

1,2, 

.  .  .  ,n 

be 

Set  x  . 

-  r  . 

=  h  . 

(A-l) 


Then  F .  (r  .  )  =  0 
i  3 

Writing  a  Taylor's  series  expansion  for  F^(r^) 


0  =  F.(r.)  =  F.(x.  -  h.) 
i  J  i  J  3 

=  F.(x.)  -  h . F ! ( x . )  +  h2  F. (x.)  -  ... 
1  v  3  3  1  J  ,1  1  J  ’ 

2  ! 


( A-2  ) 


If  hj  is  small,  that  is,  if  the  choice  of  x^  is  sufficiently 

accurate,  and  if  the  curve  F.(x.)  does  not  have  derivatives 

at  or  near  x^  so  large  that  they  interfere,  then  an 

approximate  value  of  h.  obtained  by  neglecting  all  but  the 

J 

first  two  terms  in  the  series  is 


’ 

T-<  i 

-  ' 


' 


(45) 


Then  a  new  approximate  value  for  x.  becomes 


x  ( 2  ^  =  x(])  -  F . ( x ( 1 ^ 

J  J  i  ,] 


(A-4) 


f!  (xw 

1  J 


This  procedure  is  repeated  to  obtain  a  better  value  for  xj 


(3) 


(3)  _  (2) 


X  . 

J 


or  in  general 


F  (x(2) 
i  .1 

f'(x(2) 

i  j 


(A-5) 


x(r)  =  XU-1)  _ 

3  J 


F.  (x 

i  ,1 


(r-1 ) 


F'(xlr-1) 

i  j 


The  procedure  converges  to  a  solution  for  x.  when 

J 

(r)  ( r-1 ) 

x.  '  -  x.  '  +  specified  tolerance. 

J  J 

For  calculation  purposes  an  n  dimensional  problem  becomes 

F.hx^-1))  x*4*  =  x^"1)  -  F.lx^:1)) 

1  J  J  j  J  J  i v  J-l  ’ 

or  in  matrix  form 


f’U^-1))  .  x<d>  =  f’U^"1))  x^-D  -  Ffx^-D) 


(A-6) 


' 


.  • 


. 


«*.  MUI*!  '  lo4  t>0  i  t  i  »<]<*  ±.  1  .<  ia  ^  X 


' •  •  /’< »  v,.  ^;Xilfc^ 

I 
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Nomenclature 


F  .  ( x  .  ) 
iv  j' 


F"(x . ) 
v  j' 


empl oyed 

equation  to  be  solved;  i  =  l,2,...,n; 
first  derivative  of  F.(x.) 

second  derivative  of  F.(x.) 

iv  3J 

first  approximation  to  true  solution 

t  h 

d  approximation  to  true  solution 


j  =  1,2,, 


(ii)  Derivation  of  Finite-Difference  Equations 


By  using  equations  (3-15)  to  (3-18)  inclusive  it  is  possible 
to  derive  the  general  finite-difference  equations  representing 
the  energy  balance  for  spherical  coordinates. 

For  n  =  1 , 2 , 3 , . . . , n-1 ,  the  energy  equation  is; 


(A-7) 


T  ,  -  2  T  +  T  ,  +2  ( 2  T  -.+3T  -  6T  -,+T  n 

n+1  n  n-1  v  n+1  n  n-1  n-2 


A  r‘ 


m^r 


6Ar 


-  r  Hifi(cA’cB’ • • • >T>  =  0 


For  n  =  N,  the  energy  equation  is; 


* 


' 


■ 


■ 


.('*7) 


Tn_2  +  10Tn_1  -  25Tn  +  16Ts  +  2  (3Tn_2  -  20T 


N-l 


.  ( A-8 ) 

15Tn  +  32Tg) 


5a  r 


2 


mjr 


30Ar 


-  >_  Hifi(CA,C  .  .,T)  =  0 

i 


where  mAr  =  r 


By  using  equations  (3-19)  to  (3-22)  inclusive,  we  can  derive 
the  general  finite-difference  equations  representing  the 
mass  balance  in  spherical  coordinates. 


For  n  =  1 , 2 , 3 , . . . , n-l ,  the  mass  equation; 


(a-9  ) 

C.  ,  -  2C.  +  C.  ,  +  2  ( 2C .  ,  +  3C.  -  6C .  ,  +  C.  0 

l .  n+1 _ l ,  n _ l ,  n-l  _ v  l ,  n+1  i  ,  n _ l ,  n-l _ l  .  n-2 


) 


hr 


mAr 


6  hr 


-  f.(CA,CB, - T)  =  0 


For  n  =  N,  the  mass  equation  is; 


C.  XT  0  +  IOC.  XT  n  -  25C.  XT  +  l6C .  + 

i , N-2  l , N-l  l  , N  l  ,  S 

cr *  2 

5at 


(a— io ) 


2  <5Ci.N-2  -  20C i , N— 1  "  15C1,N  +  32Ci.S> 


mAr 


30  Ar 


fi(CA>CB>  •  •  •  ’T)  =  0 


. 


- 


. 


\ 


v  ;4  '  . 


('■8) 


(iii)  Calculation  of  Overall  Effectiveness  Factor 


We  want  to  numerically  evaluate  the  integral; 


( A-l 1 ) 


n 


R 

3r^  <5  dr 

R3  O 

o  R 


where  0  is  the  rate  of  reaction  at  grid  point  n. 


Assume  we  want  to  approximate  the  function  by  passing  ‘ 
a  fourth  degree  polynomial  through  any  five  known  points 
on  the  radii.  This  requirement  enables  one  to  develop  a 
formula  which  has  the  following  form; 


N 


N 


(A-12) 


f(x)dx  = 
o  0 


W  f ( x  )  +  ...  + 

n  v  n ' 


W  f  ( x  ) 
n  v  n ' 


n=N-4 


where  W  are  the  constants  to  he  determined, 
n 


For  any  set  of  known  five  grid  points  using  a  total 
interval  of  one  with  x^  =  r^  we  obtain  the  following 
linear  matrix; 


L(l)  =  Wj  '  +  W2  +  W3  +  W4  +  W5 


=  1 


(A-13) 


L(x)  =  +  W2r2  +  W5r3  +  W,(r4  +  W  ?r  ? 


=  1 


. 


. 


. 


. 


■ 
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L  ( x2 ) 

=  "1 

2 

ri  + 

to 

to  to 

+ 

W  r2 
3  3 

+ 

T,  .  2 
W4r4 

+  W-rj:  =  1 

P  2  3 

L(x5) 

=  >i 

r3  + 

1 

W  r3 

2  2 

+ 

W  r 3 
vv3i3 

+ 

W  r 3 
vv4r4 

+  W  r3  =  1 

4 

L(x4) 

=  wi 

4 

rl  + 

W  r^ 

2  2 

+ 

W 

3  3 

+ 

W.  r? 
4  4 

+  Wcr^  =  1 

5  5  5 

For  uni  for 

m  grid 

spac ing 

r ^  to 

r 

5 

take 

on  the  values 

0,  1  , 

1  , 

3  ,  and  1  r 

espectively . 

For 

a  grid  spacing 

4 

2 

4 

consisting 

of  a 

one-half 

unit 

and 

three 

even  units 

r  to  r  take  on  the  values  0 ,  1_  ,  2  »  5.  >  and  1  * 

5  111 
respectively. 

Let  A  represent  f(r^k),  j  =  1,2,... ,5;  k  =  1,2,. ..,5 

W  represent  the  required  weighting  factors 

b  represent  RP ,  P  =  1,2,... ,5;  R  =  1 

P 

(A-14) 

Then  V  =  b  A-1 


When  W  is  obtained,  its  elements  are  multiplied  by  the 
total  interval  which  the  corresponding  five  points  represent. 
Then  since  W^m ^  and  represent  the  same  grid  point, 

they  are  combined.  The  weighting  factors  are  now  in  the 
form  in  which  they  can  be  used  for  the  problem  we  are 
about  to  solve.  Equation  (A-ll)  then  becomes; 


.  ,  I  W.  -  (c  zhll 

•* 


, 


' 


. 

. 


. 

. 
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n 


o 


3_  (  +  W^>2 

R  v) 


R 


r2  +  •••  +  V>NRN  > 


til 


(A-15) 


where  subscript  N  indicates  the  n  vector. 


(iv)  Temperature  Dependent  Rate  Constant 

The  expression 
-E 


,  mm  R  T 
ky  T  e  g 


0  <  m  <  1 


(A-16) 


summarizes  the  predictions  of  the  simpler  versions  of  the 
various  theories  for  the  temperature  dependency  of  the  rate 
constant.  Because  the  exponential  term  is  so  much  more 
temperature  -  sensitive  than  the  Tm  term,  the  variations 
of  k  caused  by  the  latter  is  effectively  masked,  this 


results  in 


k  =  k  e 
v  s 


-E 


R  T 


(A-17) 


•  For  computational  purposes,  k  must  be  eliminated,  the 

s 

reason  being  that  the  computer  cannot  accurately  calculate 

(30) 


the  difference  of  two  number  of  the  order  of  e 


(A-18) 


thus 

k 


E 


E 


k  e 
vs 


R  T 
g  s 


R  T 
g 


. 

' 


. 

. 
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B  NUMERICAL  EXAMPLE  PROBLEM 

To  show  the  applicability  of  Parameter-Perturbation  using 
Newton's  method,  we  chose  to  solve  a  non-linear  system  of 
algebraic  finite-difference  equations  found  in  detail  in 
the  book  by  Lapidus,  (8).  Lapidus  made  use  of  a 
quadratic  forcing  technique  to  obtain  solutions  to  the 
finite-difference  equations,  but  the  rate  of  convergence 
using  this  technique  is  very  slow. 

i 

Following  is  an  outline  of  the  problem. 

Given  the  finite -difference  representation  of  a  second- 
order  reaction  taking  place  in  a  fluid  flowing  through  a 
tubular  reactor  in  laminar  flow,  we  want  to  solve  for  the 
concentration  of  the  reactant  as  a  function  of  the 
dimensionless  distance. 

The  boundary  conditions  are : 


1  = 

f  -  1  ( 

0  Peh 

L  - 

f  )  at  Z 

o  ' 

=  0 

fN 

=  fN+l 

at  Z 

=  1 

where 

• 

• 

f . 

l 

represents 
i  =  1,2,.. 

the 
•  ,N 

concentration 

of  the 

reactant , 

h 

represents 

the 

dimensionless 

dis tanc 

e  segment 

N 

represents 

the 

total  number  o 

f  grid 

points 

Pe 

represents 

the 

Pec  let  number 

Z 

represents 

the 

dimensional  distance 

' 


. 

■ 
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The  specified  input  parameters  are: 
h  =  .04 
N  =  24. 

Pe  =  1.0 
R  =2.0 

where  R  is  the  Reynolds  number 

A  system  of  twenty-four  algebraic  finite-difference 
equations  having  the  following  form  can  be  generated. 

(B-2) 

2  _  — f  +  RPeh2  f2  -  (l-Peh/2)  fp  -  -1  +  Pe.h/2  =  o 
1+Peh  1  1  +  Peh 

-  ( 1+Peh/ 2 )  fj[  +  2.0  f2  +  RPeh2  f2  -  (l-Peh/2)  f^  =  0 


-  (l+Peh/2)  fn_1  +  2.0  f  +  RPeh2  f2  -  (l-Peh/2)  f  =0 


-  (l+Peh/2)  f25  +  ( 1+Peh/ 2 )  fg4  +  RPeh2  f24  =  0 

This  system  of  twenty-four  non-linear  equations  were  solved 
on  the  I.B.M.  7040  computer  using  Parameter-Perturbation 
with  Newton's  method.  Trial  vectors  used  were  F^°^=  (l.O) 
and  F^°'=  (0.0)  .  For  F^°^=  (l.O)  Newton’s  method  by 

itself  was  satisfactory,  for  F^°^=  (0.0)  one  set  of  ’perturbed' 
equations  had  to  be  solved  before  the  true  solution  to  the 
system  was  obtained.  A  comparison  of  results  obtained  by 


. 


t 


. 


. 


■. 


-(■53) 


Parameter-Perturbation  versus  results  obtained  by  Lapidus 
is  presented  in  Table  4,  page  53*  Complete  results  are 
found  in  Table  5>  page  54  * 

TABLE  4  Comparison  of  Results,  Example  Problem 


Grid  Parameter-Perturbation 


Point 

Results 

5 

.577403 

10 

.528069 

15 

.492703 

20 

.471153 

Lapidus  'Results 

%  Difference 

.577406 

t 

.0005 

. 528072 

.0006 

.492705 

.0004 

.471155 

.0004 

as 


where  *}c  difference  is  calculated 

Lapidus '  Results  -  Parameter-Perturbation  Results 

Parame ter -Per turbat ion  Results 


X 


100 


■ 


- 


■ 


TABLE  5 


Complete  Results,  Example  Problem 


Trial  Number  of  Sets  of 

Vector  Perturbed  Equations 

Required 


F  ^  °  ^  =  1.00 


0 


F  ^  °  ^  =  0.00 


1 


Solution  to 
the  Set  of 
Fini te-Ui f f erence 
Equations 

.  627874 
.614265 
.601332 
.589052 
.577403 
. 566368 

.555930 
.546075 
.536791 
. 528069 
.519902 
.  5H228 
.505212 
.498685 
.498685 
.492703 
.487270 
.482390 
.478070 
.474321 
.471153 
.468581 
.  466620 
.465290 
.465290 

Same  as  those  for 


'(54) 


Object  Time 
IBM  7040 


4  seconds 


seconds 


1.000 


i 

. 

m 


■ 


■ 

■■■'■*  X  *  '  '  '  I 


(55) 


From  Table  4,  page53  and  Table  5,  page  54  it  is  noted  that, 

(1)  Newton's  method  increases  the  rate  at  which  convergence 
is  obtained  as  compared  to  the  801  iterations  required 
by  the  quadratic  forcing  technique. 

(2)  Parameter-Perturbation  enables  one  to  use  any  physically 
possible  initial  trial  vector. 

(3)  Calculation  time  is  small,  4  to  7  seconds  on  an 

IBM  7040  computer,  compared  to  Lapidus  method  which, 
took  40  seconds  per  100  iterations  on  an  IBM  701 
computer . 


► 
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C  GENERAL  FLOW  DIAGRAM 


Start 


Define  N,  N1 ,  N2,...,N9,  function  of  the 
total  number  of  finite-difference  equations 


i 

> 

J  Define  system  input  parameters;  DALP,  DELH 
J  DELP,  DEX,  DKV,  DKV1,  PX,  R,  TS,  W 


Define  "True"  equations  P.^  =  f(system 
parameters),  i  =  1,...,N,  k  =  l,..,Z+2 


Define  ’Trial'  vectors,  ,  i  =  1,.,N 


;  Define  'Derived'  equations  Q.,  ,  where  Q . , 
P.,  except  when  i=k,  i=l^...,N 


_ t _ _ _ _ 

Set  integer  counters  equal  to  1 


Is  Dljg  >  1.0  ?  No 


Yes 

< 


Set  D1i2  =  i.o 


A 


' 


. 


. 


■  ■  •  ;  «' ;  f 


. 


. 


- 
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!  40 


Calculate  Ql.,,  where  Q1  Is  defi 

~  «hk  =  <4k  +  (pik  -  <iik)ki>iI2/v 


ned  as 


100 


41 


Evaluate  the  partial  derivatives, 
^ik  =  ( ^ ( Per turbed  equations ) ^ ) /x. 


i 


Evaluate  the  residuals  ,  i  =  1,...,N 


Calculate  and  store  the  absolute  sum  of 
the  residuals 


Calculate  R, 

I 

H.S.  of  matrix  A 

r 

Increase  11  by  1 

i  _  . 

.  _ _ _  i 

• 

L  _ 

Solve  for  X.n+^ 
i 

i  1 

- . - _ _ i 

Evaluate  the  residuals  ,  i  =  1,...,N 


Calculate  and  store  the  residuals 


' 


f 


. 


. 

■ 

t 

' 

■ 


■ 

■ 


¥ 


(58) 


3 

L..  i 


Is  C4  *  C4  ? 


No 


*  67 


Yes 

Increase  12,  13  by  1 


Is  13  >  12  ? 


Yes  >.888 


"'■r 

No 


i 


Set  X.  =  F. 

l  l 


;  888  j 
■-L . — 1 


>  Increase  15  by  1 


1 


4 


* 


i  '  '  - 


■ 


■ 


■ 
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h 


t 

| 


I 


Go  to  100 

1 _  _ _ _ 


I 


I  67  I - H  Define  maximum  allowable  residual,  V 


1 

Is  GiO  ? 


Yes  ^  99 

I 


No 


» 


* 


(.60) 


99 


Store  current  value  of  D1j9,  set  B  =  Dl^ 


is  W  ? 


Yes  _  _>*  150 

l _ i 


No 


Set  F±  =  Xi 


Increase  12  "by  1 


!  Set  C6 


V  __ 

=  ABs  (D1 


12-2 


D1 12-1 ^ 


i 


I  100 


150 


Calculate  rate  of  reaction  at  surface  and 
all  interior  points  R.  j  =  0,1, ...,Z 

J 


5 


1 


. 


. 


V 


. 
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I 

Calculate  local  effectiveness  factor, 

EFF  .  =  R ./R 
J  3  0 

j 

I 

Evaluate  overall  effectiveness  factor  EFF 

_ 

1 

Evaluate  GAMMA,  BETA,  THMOD 


1 

Write  : 

(1) 

GAMMA , 

,  BETA,  THMOD 

(2) 

EFF  .  , 

J 

j  =  0, 1 , . . . ,Z 

(3) 

EFF 

(M 

V  1 

=  1 , 2, . . . , 12 

(5) 

xi’  1 

=  1,2,  ...  ,N 

(6) 

PX 

• 

(7) 

cxs 

(8) 

DELR 

(9) 

DELH 

(10) 

TS 

(11) 

Gi’  1 

=  1 , 2 ,  .  .  .  ,  N 

< 


L 


1 


■ 


* 


■  •  •)  ,  ; 

' 


I' 


(.62) 


i 

200  -c  No 


Have  all  input  parameters  been  used? 


A 

Yes 

* 


Stop 


. 


1  l. 


V 


.  :•  V 


•  i 
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Program 

Au 

B 

B1 

BETA 

CXS 

C4(ll) 

DALP 

DELE 

DELR 

DEX 

I)KV 

DKV1 
El (12) 
EAC 


Nomenclature 

Partial  derivatives,  I  = 

Value  of  increment  Dl(l2) 
of  convergence 

Present  value  of  15,  in  f 

Dimensionless  parameter, 
CXS,  DEX,  DALP,  TS) 

Surface  concentration  of 

Sum  of  residuals 

Thermal  conductivity  of  c 

Heat  of  reaction 

Incremental  radii,  square 

Effective  diffusion  const 

Rate  constant  at  surface 

Effective  adsorption  cons 
Incremental  value,  O^Dl 
Activation  energy 


1 ,  .  .  .  ,  N  ;  J  =  1 ,  .  .  .  ,  N 

,  most  recent  point 

ixed  mode 

function  of  (DELH, 

i 

X  component 

atalyst  particle 

d 

ant,  component  X 

conditions 

tant 
(l2)£  1 


. 


i  '  ■ 


■ 


■ 


. 


(6'i ) 


EFF 

Overall  effectiveness  factor 

EFF(l) 

Local  effectiveness  factor(s),  I  =1,  2,  Z 

EK 

Defined  as  EAC/(RG  TS) 

Fi 

Most  recent  correct  solution  to  'perturbed' 
equat  i  oils 

GAMMA 

Same  as  EK 

W 

>"3 

M 

Integer  variables  1 

IK, 12, 13, 15 

Integer  counter,  used  for  convergence 
procedure 

LI,  L2 ,  •  •  « 

Integer  variables,  used  for  parameter 
variations 

N 

Number  of  f inite-dif f erence  equations 

PIK 

Parameter  of  finite-difference  equations 
representing  the  physical  problem,  1=1,2, 
...,  N ;  K  =  1,  2,  Z+2 

PC 

Proportionality  constant 

PT 

Total  pressure  of  system 

PX 

Partial  pressure  of  X  component 

qik 

Parameters  of  'derived'  equations,  1=1,  2, 
...,N;  K=l,  2,  ...,  Z+2 

. 

■ 


t'li 


. 


' 


■ 
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Qi 


IK 


RG 

R 

THMOD 

TS 

V 

¥ 


Z 


Parameters  of  'perturbed'  equations,  defined 
as  Qlr,  K  =  Qik  +  (PJK  -  Qik)  D1(I2)/W, 

1=1,  2,  .  .  . ,  N;  K  =  1,  2,  . . . ,  Z+2 

Gas  constant,  either  1.087  or  82.1 
Radii  of  spherical  catalyst 

Modified  Thiele  Modulus 


Surface  temperature 

Specified  maximum  residual  size 

Maximum  incremental  value 


Approximate  solution  to  present  set  of 
'perturbed'  equations 

Number  of  grid  points  used  in  finite- 
difference  representation 


. 
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D.  TABULATED  RESULTS 

(i)  First  Order  System 


Type  of  Computer  IBM  7040 

Number  of  Sets  of  Solutions  Obtained  22 

Number  of  Grid  Points  Employed,  Spherical  23 

Object  Time  9  min. 

Fortran  Time  3  min. 


45 

41 


Type  of  Computer 

Number  of  Sets  of  Solutions  Obtained 
Number  of  Grid  Points  Employed,  Spherical 
Object  Time 
Fortran  Time 


IBM  7040 
24 
19 
10 
3 


min^  1 
min.  26 


(ii)  Langmuir-Hinshelwood  System 


Type  of  Computer 

Number  of  Sets  of  Solutions  Obtained 
Number  of  Grid  Points  Employed,  Spherical 
Object  Time 
Fortran  Time 


IBM  7040 
12 
23 
8 
3 


min . 
min . 


52 

3'* 


Type  of  Computer 

Number  of  Sets  of  Solutions  Obtained 
Number  of  Grid  Points  Employed,  Spherical 
Object  Time 
Fortran  Time 


IBM 


7040 

15 

19 

5 

3 


min . 
min . 


34 

6 


Type  of  Computer 

Number  of  Sets  of  Solutions  Obtained 
Number  of  Grid  Points  Employed,  Slab 
Object  Time 
Fortran  Time 


IBM  7040 
15 
19 

6 

3 


min . 
min . 


18 

10 


sec . 
sec  . 


sec 

sec 


sec 

sec 


sec 

sec 


sec 

sec 


'  I  l  (  •  ) 
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■ 
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(ill)  Example  of  Computer  Output 


System  Parameters 


Gamma 

30.0  , 

,  Beta 

0.00  , 

Th .Mod . 

1.79 

,  P 

xs 

10.0 

C 

.003  . 

,  Delr 

.0019  , 

Delh 

0.00 

,  T 

402.6 

xs 

s 

Dex 

.005  ; 

„  Dalp 

.0005  , 

Dkv 

.01 

Overall  Effectiveness  Factor  .88901 
Number  of  Increments  -1.00 


Position 

Local  Effectiveness 

Concentration 

Teinperatur 

r 

Factor 

C 

'  T 

R 

C 

T 

s 

s 

1/46 

.7309 

.73097 

1.00000 

3/46 

.7319 

.73189 

1.00000 

5/46 

.7337 

.73373 

1.00000 

7/46 

.7365 

.73651 

1.00000 

9/46 

.7402 

.74021 

1.00000 

11/46 

.7448 

.74486 

1.00000 

13/46 

.7504 

.75046 

1.00000 

15/46 

.7570 

.75703 

1.00000 

17/46 

.7646 

.76458 

1.00000 

19/46 

.7731 

.77313 

1.00000 

21/46 

.7827 

.78269 

1.00000 

23/46 

.7933 

.79329 

1.00000 

25/46 

.8049 

.80496 

1.00000 

27/46 

.8177 

.81771 

1.00000 

29/46 

.8316 

.83159 

1.00000 

31/46 

.8466 

.84662 

1.00000 

33/46 

.8628 

.86283 

1.00000 

35/46  . 

.8803 

.88027 

1.00000 

37/46 

.8990 

.89897 

1.00000 

39/46 

.9190 

.91898 

1.00000 

41/46 

.9403 

.94034 

1.00000 

43/46 

.9631 

.96311 

1.00000 

45/46 

.9873 

.98733 

1.00000 

Residuals 

<  .00001 

-  1 


!  .  -  ' 


. 


. 
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)  Summary  of  Results 

ctiori  Mechanism 

R 

J 

A. 

n 

First  Order 

0.00 

30.0 

1.778 

.889 

(Spherical  - 

0.00 

30.0 

4.467 

.745 

Coordinates ) 

0.00 

30.0 

6.537 

.654 

0.00 

30.0 

8.263 

.590 

0.00 

30.0 

9.769 

.543 

0.00 

30.0 

11.120 

.505 

0.00 

30.0 

12.356 

.475 

0.00 

30.0 

13.501 

.450 

0.00 

30.0 

14.574 

.429 

0.00 

30.0 

15.587 

.410 

-.20 

30.0 

1.275 

00 

VO 

« 

-.20 

30.0 

2.736 

.456 

-.20 

30.0 

3.794 

.380 

-.20 

30.0 

4.669 

.333 

-.20 

30.0 

5.432 

.302 

-.20 

30.0 

6.119 

.278 

-.20 

30.0 

6.746 

.259 

-.20 

30.0 

7.329 

.244 

-.20 

30.0 

7.875 

.232 

-.20 

30.0 

8.390 

.  221 

+  .05 

30.0 

2. 108 

1.054 

+  .05 

30.0 

6.006 

1.001 

. 


’ 


«  I 


' 
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Reaction  Mechanism  P. 

A 

Langmuir-  r- . 

Hinshe lwood 

**"  • 

(Spherical  - 
Coordinates) 


(Linear  - 
Coordinates  ) 


$  f) 


7.79s 

.135 

4.073 

.  212 

2.968 

.251 

2.396 

.291 

2.036 

.318 

3.020 

.0524 

1.750 

.0910 

1.330 

.  116 

1.110 

.135 

.969 

.151 

IC 

s  o 

96  4 

964 

964 

964 

964 

964 

964 

964 

964 

964 


. 


(v)  Fortran  Source  List 


Following  is  a  summary  of  a  general  program  utilizing 
Parameter-Perturbation  to  solve  the  rate  expression 

r.  =  k . P . /( 1+K . P . ) 

Program  1  is  for  spherical  coordinates,  Program  2  is 
for  a  slab.  All  input  parameters  must  be  specified 
as  well  as  the  number  of  grid  points  to  be  employed. 

Since  Programs  1  and  2  are  quite  similar,  Program  2 
consists  of  only  a  list  of  cards  which  are  different. 
Unless  stated  otherwise  assume  the  cards  are  interchange 


able  . 


•  ' 
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MB 

1 

2 

ry 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 
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STATEMENT 

NUMBER 


1 


STATEMENT 


GRID  =  15  OR  19  OR  23  OR  .... 

N4  =  GRID 
N1  =  2  *  N4 
N2  =  N4  +  2 
N3  =  N4  +  1 
N5  =  N4  -  1 
N6  =  N1  -  1 

N7  =  N4  +  3 
N8  =  N1  +  1 
N9  =  N4  -  2 
NN1  =  N1  -  2 
FORMAT  (IX, 514.8) 

READ  (5,1 )TS,DELH,EAC,RG,PX 
READ  ( 5 , 1 ) PC , DALP , DEX , DKV , DKV 1 
CXS  =  PX/(RG  *  TS) 

EK  =  EAC/( 1.987  *  TS) 

DK0  =  DKV 

Z  =  PC  *  DEX/(RG  *  TS) 

DN1  =  N1 
DN6  =  N6 
DR  =  2./DN1 
DELR  =  DR  *  *  2 

DIMENSION  P ( N 1 , N 2 ) ,  X(Nl),  F(Nl) 

P(l,l)  =  Z* ( -1 ./DELR+2 . *DN1* ( -3-/(6. *DR ) ) ) 
P(l,2)  =  -P(l,l) 

P ( 2 , 1 )  =  Z*(l. /DELR+2. *DNl/3. *(-5. /(6. *DR) )) 
P(2,2)  =  Z*(-2. /DELR+2. *DNl/3. *(3. /(6-*DR) )) 
P(2,3)  =  Z*( 1. /DELR+2. *DNl/3.*(2./( 6. *DR) ) ) 
D0  508  I  =  3,  N5 


v  • 


■ 


• 
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CARD  STATEMENT  STATEMENT 


NUMBER 

NUMBER 

30 

is 

i— 1 

m 

i— i 

31 

Z2  =  2 . *Z 1  -  1. 

32 

P  ( 1 , 1-2  )  =  Z*(2.*DNl/Z2*(  1  .'/(6. *DR)  )  ) 

33 

P(l,I-l)  =  Z8(l./DELR+2.*DNl/Z2*(-6./(6.*DR) ) ) 

34 

P(I,I)  =  Z*(-2./DELR+2.*DNl/Z2*(3./(6.*DR) ) ) 

35 

P(I,I+1)  =  Z*(l./DELR+2.*DNl/Z2*(2./(6.*DR) ) ) 

36 

508 

CONTINUE ' 

37 

P(N4,N9)  =  Z*(-l./(5.*DELR)+2.*DNl/DN6*(3. 

/( 30 . *DR )  )  ) 

38 

P(N4,N5)  =  Z*(l0./(5.*DELR)+2.*DNl/DN6*(-20 
./( 30 . *DR) ) ) 

39 

P(N4,N4)  =  Z* (— 25  » / (5. *DELR ) +2 . *DNl/DN6* ( - 
15-/(30.*DR)  )) 

40 

B0  557  I  =  N3,  N1 

41 

B0  557  J=  1,  N4 

42 

LI  =  I  -  N4 

43 

557 

P(I,J)  =  P(L1 , j)*ts*dalp*rg/(pc*dex) 

44 

P(N4 ,N2 )  =  Z*(l6./(5.*DELR)+2.*DNl/DN6* 

( 32 ./ ( 30 . *DELR )  )  )*PX 

45 

P(N1,N2)  =  TS*DALP*(l6./(5.*DELR)+2.*DNl/ 
DN6*(32./(50.*DR)  )  ) 

46 

P(1,N3)  =  -  PC*DKV 

47 

D0  701  K  =  2,  N4 

48 

701 

P(K,N3)  =  P(1,N3) 

49 

P(N3,N3)  =  DKO  *  DELH 

50 

J)0  703  K  =  N 2 ,  N1 

51 

703 

P(K,N3)  =  P(N3,  N3) 

52 

D0  197  I  =  1.  N4 

53 

197 

F(l)  =  PX 

54 

D0  188  I  =  N3,  N1 

55 

188 

F(l)  =  1.00000 

56 

D0  157  I  =  1,  N1 

i  e  * 


. , 
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CARD 

.STATEMENT 

STATEMENT 

NUMBER 

NUMBER 

57 

157 

X  (i)  =  F(I) 

58 

DIMENSION  Q(N1,N2),  Ql(Nl,N2),  Dl(50) 

59 

D0  107  1=1,  N1 

60 

D0  107  J  =  1,  N2 

6l 

107 

Q  ( I ?  J )  =  P(I,J) 

62 

Q  ( 1 , 1  )  =  (Q ( 1 ?  2 ) *X( 2 ) +Q ( 1 , N3 ) *X( 1 )/ 
(l.+DKVl*X(i)))/x(i) 

63 

Q  (  2 , 2  )  =  -(Q(2,l)*X(l)+Ql(2,3)*X(3) 

+Q  ( 2 ,N3  ) *X( 2 )■/(}..  +DKV1*X(2  ) )  )/x( 2 ) 

64 

B0  110  I  =  3,  N5 

65 

Q(I,I)  =  -  (Q(l.I-2)*X(l-2)+Q(l,I-l) 
*X(l-L)+Q( 1,1+1 )*X( 1+1 )+Q (i ,N3)*X(l )/ 

( 1 . +DKV1*X( I  ) ) )/x( I ) 

66 

Q(N4,N4)  =  ~(Q(N4,N4-2)*X(N4-2)+Q(N4, 

N4-1 )*X(N4-1 )+Q (N4 ,N2 ) +Q (N4 ,N3 ) *X(N4 )/ 

( 1 . +DKV1*X(N4 ) ) )/x(N4 ) 

67 

Q(N3,1)  =  -(Q(N3,2)+Q(N3,N3)*X(1)/ 
(l+DKVl*X(l)) ) 

68 

Q(N2,2)  =  -(Q(N2,1)=Q(N2,3)+Q(N2,N3)* 

X( 2 )/ ( 1 . +DKV1*X( 2 ) ) ) 

69 

D0  112  I  =  N7 ,  N6 

70 

L2  =  I  -  N2 

71 

L3  =  I  -  N3 

72 

• 

L4  =  I  -  N5 

73  • 

L5  =  I  -  N4 

74 

112 

Q(I,L5)  =  -(Q(I,L2)+Q(I,L3)+Q(I,L4)+ 
Q(I,N3)*X(L5)/(1.+DKV1*X(L5)) ) 

75 

Q(N1,N4)  =  ~(Q(Nl,N4-2)+Q(Nl,N4-l)+ 

Q (Nl ,N4+2 ) +Q (N1 ,N3 )*X(N4 )/( 1 . +DKV1*X(N4 ) ) ) 

76 

12  =  1 

77 

13  =  1 

78 

• 

15  =  1 

79 

D1  ( 12 )  =  1.00000 

. 


- 


' 


' 


NUMB 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 


STATEMENT 

NUMBER 

100 


999 

998 

C 


30 

40 

C 


850 
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STATEMENT 


CONTINUE 

II  =  1 

IF(Dl(l2).GT.W) G0  T0  999 
G 0  T0  998 
D1(I2)  =  1.00000 
CONTINUE 

EVALUATING  THE  PRESENT  DERIVED  CONSTANTS 
B0  30  I  =  1,N1 
D0  30  K  =  1,N2 

Qi(i,K)  -  Q(i>k)+(p(i,k)~q(i,k)) *Di ( l2 ) 

CONTINUE 

EVALUATING  THE  PARTIAL  DERIVATIVES 
DIMENSION  A(N1,N8) 

A ( 1 , 2  )  =  Ql(l,2)/X(N2) 

A ( 1 , N2 )  =  -Ql(l ,2)*X(2)/X(N2)**2 
A ( 2 , 1 )  =  Q1,(2,1)/X(N3) 

A ( 2 , 3 )  =  Ql(2,3)/X(N3  +  2) 

A(2,N3)  =  -Q1(2,1)*X(1)/X(N3)**2 

A ( 2 , N3+2 )  =  -Ql(2,3)*X(3)/X(x(N3+2)**2 

D$  850  I  =  3,  N5 

L2  =  I  +  N9 

L3  =  I  +  N5 

L4  =  I  +  N3 

A ( I , 1-2 )  =  Ql(l,I-2)/x(L2) 

A ( I , 1-1 )  =  Q1(I,I-1)/X(L3) 

A ( I , L3 )  =  -Q1(I,I-1)*X(I-1)/X(L3)**2 
A ( I , L2 )  =  -Ql(l,I-2)*X(l-2)/X(L2)**2 
A(l,I  +  l)  =  Q1(I,I  +  1)/X.(L4) 

A ( I , L4 )  -  -Ql(l,I+l) *X( 1+1 )/x( L4 ) **2 
A (N4 , N4-2 )  =  Ql(N4,N4-2)/x(Nl-2) 


1 

' 
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CARD 

NUMBER 
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110 

A (N4 , Nl-2 )  =  Q 1 (N4 , N4-2 ) *X (N4-2 )/x(Nl- 
**2 

111 

A (N4 , N4-1 )  =  Q1(N4?N4-1)/x(N1~1) 

112 

A(N4,Nl-l)  =  -Q 1 (N4 , N4-1 ) *X(N4-1 ) / 
X(N1-1)**2 

113 

V0  668  I  =  1,  N4 

114 

L2  =  I  +  N4 

113 

A  ( I ,  I )  =  Q1(I.I)/X(L2)+Q1(I,N3)*( 
(1.+DKV1*X(i))*EXP(EK-EK/X(L2))- 
DKV1*EXP ( EK-EK/X ( L2 ) )*X( I ) )/ 
(l.+DKVl*X(l) )**2) 

116 

668 

A ( I , L2 )  =  Q 1 ( I ?  N3 ) * ( ( X ( I ) *EXP ( EK-EK/ 
X(L2) )*EK/x(L2)**2)/(1.+DKV1*X(I)) )- 
Q1(I,I)*X(I)/X(L2)**2 

117 

B0  669  I  =  N3,  N1 

118 

L2  =  I  -  N4 

119 

A  ( I ,  L2  )  =Ql(l;,N3)*(  (  (l.+DKVl*X(L2)  )* 
EXP ( EK  -EK/X ( I ) ) -DKV 1 *EXP ( EK-EK/X ( I ) ) 
*X(L2) )/(l.+DKVl*X(L2) )**2) 

120 

669 

•  A ( I , I )  =  Ql(l ,L2)+Ql(l , N3 ) * ( ( X ( L2 ) * 

EXP ( EK-EK/X ( I ) ) *EK/x ( I ) **  2 ) / ( 1 . +DKV 1 * 
X(L2) )  ) 

121 

A(N3,N2)  =  Q1(N3,2) 

122 

A(N2,N3)  =  Q 1 ( N  2 , 1 ) 

123 

A(N2,N7)  =  Q1(N2,3) 

124 

D/  779  I  =  N7,  N6 

125 

L2  =  I  -  N2 

126 

L3  =  I  -  N3 

127 

L4  =  I  -  N5 

128 

A ( I ? 1-2 )  =  Q 1 ( I , L2 ) 

129 

A ( I , 1-1 )  =  Q1(I,L3) 

130 

779 

A ( I , 1+1 )  =  Q 1 ( I , L4 ) 

131 

A(Nl,Nl-2)  =  Ql(Nl,N4-2) 

132 

A(N1,N1-1)  =  Ql(Nl,N4-l) 

■ 


. 

> 


■ 

e* i  .*  *  i 
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133 

134 

135 


136 

137 

138 

139 

140 

141 


142 


143 

144 

145 

146 

147 

148 

149 

150 
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C  .EVALUATING  THE  RESIDUALS 

G ( 1 )  =  Q1(1,1)*X(1)/X(N3)+Q1(1,2)* 
X(2)/X(N2)+Q1 ( 1 ,N3)*EXP(EK-EK/X(N3)) 
*X ( 1 )/( 1 . +DKV1*X( 1 )) 

G ( 2 )  =  Q1  (  2 , 1  )*  X( 1 )/x(N3 ) +Q1 ( 2 , 2 ) * 
X( 2 )/x(N2 ) +Q1 ( 2 , 3 ) *X( 3 )/x(N7 ) +Q1 

( 2 , N3 ) *EXP ( EK-EK/X (N2 ) ) *X ( 2 ) / ( 1 . + 
DKV1*X( 2 ) ) 


D0 

885 

I 

— 

L2 

=  I 

+ 

N9 

L3 

=  I 

+ 

N5 

L4 

=  I 

+ 

N4 

L5 

=  I 

+ 

N3 

885  G ( I )  =  Q1(I, I-2)*X(l-2)/x(L2)+Ql 

(I,I-1)*X(I-1)/X(L3)+Q1(I,I)*X(I) 
/X(L4)+Q1 (i, I+1)*X(I+1)/x(L5)+Q1 
( I , N3 ) *EXP ( EK-EK/X ( L4 ) ) *X( I )/ ( 1 . + 
DKV1*X ( I ) ) 

G ( N 4 )  =  Ql(N4,N4-2)*X(N4-2)/x(Nl-2) 

•  +Q1 (N4 ,N4-1 )*X(N4-1 )/X(Nl-l)+Ql(N4. 

N4)*X(N4)/X(Nl)+Ql(N4,N2)+Ql(N4,N3) 
*EXP ( EK-EK/X (N1 ) ) *X ( N4 ) / ( 1  * +DKV 1 * 

X (N4 ) ) 

G(N3)  =  Ql(N3,l)*X(N3)+Ql(N3,2)*X(N2) 
+Q 1 (N3 , N3 ) *EXP ( EK-EK/X ( N3 ) ) *X( 1 ) / 
(l.+DKVl*X(l) ) 

G(N2)  =  Ql(N2,l)*X(N3)+Ql(N2,2)*X(N2) 
+Q1(N2, 3)*X(N7)+Q1 (N2,N3)*EXP(EK-EK/ 
X(N2  )  )*X(2)/ (l.+DKVl*X(2)  ) 

B0  889  I  =  N7,  N6 


M 

II 

CM 

-  N2 

M 

II 

in 

-  N3 

L4  =  I 

-  N4 

1— 1 

11 

in 

-  N5 

889  G ( I )  =  Ql(i  L2)*X(l-2)+Ql  I,L3)*X(I-1) 

+Q1 ( I . L4 ) *X( I ) +Q 1 ( I , L5 ) *X( 1+ 1 ) +Q 1 
(I,N3)*EXP(Ek-EK/x(I)*X(L4)/(1.+DKV1* 
X(L4 ) ) 


' 

. 


* 


■ 


. 
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151 

G(N1)  =  Ql(Nl,N4-2)*X(Nl-2)+Ql(Nl, 
N4-1 )*X(N1-1)+Q1 (N1 ,N4)*X(N1 )+Ql 
(N1,N2)+Q1(N1,N3)*EXP(EK-EK/X(N1) ) 
*X(N4)/ (1 ,+DKVl*X(N4) ) 

152 

C 

SUMMING  ABSOLUTE  VALUES 

153 

DIMENSION  C 4 ( 5 0 ) ,  R(N4),  EFE(N4) 

154 

G1  =  0.0 

155 

D0  600  I  =  1,  N1 

156 

G2  =  ABS  (G(l) ) 

157 

600 

G1  =  G2  +  G1 

15S 

C4 ( 11 )  =  G1 

159 

11  =  11  +  1 

160 

C 

EVALUATING  THE  RIGHT  HAND  SIDE  OF 
EQUATIONS 

161 

A(i,N8)  =  -G ( 1 ) +A (l.l) *X ( 1 ) +A (l,2)* 
X( 2 ) +A ( 1 , N2 ) *X (N2 ) +A ( 1 , N3 ) *X (N3 ) 

162 

A ( 2 , N8 )  =  -G ( 2 ) +A ( 2 , 1 ) *X( 1 )  +A (2,2) 
X ( 2 ) +A ( 2 , N2 ) *X ( N2 ) +A ( 2 , 3 ) *X ( 3 ) +A ( 2 , 

N3)*X(N3)+A(2»N2)*X(N2) 

163 

B0  883  I  =  3,  N5 

164 

L2  =  I  +  N9 

165 

L3  =  I  +  N5 

166 

L4  =  I  +  N4 

167 

L5  =  I  +  N3 

168 

883 

A(I,N8)  =  -G(I)+A(I, I-2)+X(l-2)+A 
(I,I-1)*X(I-1)+A(I.I+1)*X(I+1)+ 

A ( I , I ) *X( I ) +A ( I , L2 ) *X ( L2 ) +A ( I , L3 ) 
*X( L3 ) +A ( I , L4 ) *X( L4 ) +A ( I , L5 ) *X( L5 ) 

169 

A  f N4 , N8 )  =  -G(N4)+A(N4,N9)*X(N9)  + 
A(N4,N5)*X(N5)+A(N4,Nl-2)*X(Nl-2) 
+A(N4,N1-1)*X(N1-1)+A(N4,N4)*X(N4) 
+A(N4,N1)*X(N1) 

170 

A(N3,N8)  =  -G(N3)+A(N3,1)*X(1)+A( 

N3 , N3 ) *X(N3 ) +A (N3 , N2 ) *X(N2 ) 

171 

A (N2 , NS )  =  -G(N2)+A(N2, 2)*X(2)+A 
(N2 , N3 ) *Xf N3 ) +A (N2 , N2 ) *X(N2 ) +A 
(N2,N7 )*X(N7 ) 

- 

■ 
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173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 
199 
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D0  767  I  =  N7  j  N6 
L2  =  I  -  N4 

767  A  ( I ,  N8  )  =  -G(l)+A(l,L2)*X(L2)+ 

A(I, I-2)*X(l-2)+A(l , 1-1 )*X(I~1) 

+a(i,i)*x(i)+a(i>i+i)*x(i+i) 

A(N1,N8)  =  -GCNlJ+AfNl^^X^) 
+A(N1,NN1)*X(NN1)+A(N1,N6)*X(N6) 

+a(ni,ni)*x(ni) 

C  FINDING  THE  NEW  IMPROVED  ROOTS 

N  =  N1 

19  M  =  N+l 

K1  =  N 

13  C  =  A(K1,K1) 

D0  91=1)  K1 

9  A (K1 , I )  =  A(K1,I)/C 

A(K1,M)  =  A(K1,M/C 
I  =  K1  -  1 

11  D  =  A(I,K1) 

DJZf  10  J  =  1,  K1 

10  A(I,J)  =  A(I,J)  -  D*A(K1,J) 

A(l,M)  =  A(I,M)  -  D*A(K1,M) 

1  =  1-1 
L  =  0 

IF(l.EQ.L)  G$  T0  12 
G0  T0  11 

12  K1  =  K1  -  1 
L  =  1 

IF(Kl.EQ.L)  G0  TO  14 
G$  T0  13 

14  E  =  A(K1,  Kl) 

A(K1,M)  =  A(K1,M)/E 


SUM 


0 


■ 


. 

. 
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200 

J)0  1 6  I  =  2 ,  N 

201 

t— 1 

ii 

w 

1 

i— 1 

202 

D0  15  J  =  1,  K1 

203 

15 

SUM  =  SUM  +  A(I,J)*A(J,M) 

204 

A(I,M)  =  A(I,M)  -  SUM 

205 

16 

SUM  =  0 

206 

D0  18  I  =  1,N 

207 

18 

X(  I )  =  A ( I , M) 

2208 

C 

EVALUATING  THE  RESIDUALS  WITH  NEW 

ROOTS 

209 

6  ( 1 )  =Q  1  ( 1 ,  i  )  *X(  i,)/x(N3  )  +Q 1  ( 1 , 2  )  *X(  2  ) 

/X (N2 ) +Q 1 ( 1 , N3 ) *EXP ( EK-EK/X(N3 ) ) * 

X( 1 )/( 1 „ +DKV 1*X( 1 ) ) 

210 

G(2)=Ql(2, l)*X( 1)/X(N3)+Ql(2 , 2)*X(2)/ 
X(N2 ) +Q 1 ( 2 , 3 ) *X( 3 )/x(N7 ) +Q 1 ( 2 , N3 ) * 

EXP  (EK-EK/x(N2  )  ) *X(  2  )/( 1 .  +DICV  1*X(  2  )  ) 

211 

DO  886  1=3, N5 

212 

L2=I+N9 

213 

L3=I+N5 

214 

L4=I+N4 

215 

L5=I+N3 

216 

886 

G ( I ) =Q 1 ( I , 1-2 ) *X( 1-2 )/x( L2 ) +Q 1 ( I , 1-1 ) 
*X(I-1)/X(L3)+Q1(I,I)*X(I)/X(L4)+ 
Q1(I,I+1)*X(I+1)/X(L5)+Q1(I,N3)*EXP( 
EK-EK/x( L4 ) ) *X( I )/( 1 . +DKV1*X( I ) ) 

217 

G (N4 ) =Q 1 (N4 , N4-2 ) *X(N4-2 )/x(N1-2 ) +Q 1 ( 
N4,N4-1)*X(N4-1)/X(N1-1)+Q1(N4,N4)*X(N4) 

/X(N1 )  +Q1  (N4  ,N2  )  +Q 1  (N4  ,  N3  )  *EXP  (EK-EIC/ 
X(N4)/(l.  +DK\rl*X(N4  )  ) 

' 


'  ' 

■  •  •  /  i 
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218 

G(N3)=Ql(N3 , l)*X(N3)+Ql(N3,2)*X(N2) 

+Q 1 (N3 , N3 ) *EXP ( EK-EK/X (N3 ) ) *X ( i )/ 

( 1 . +DKV 1*X( 1 ) ) 

219 

G(N2 )=Q1 (N2 , l)*X(N3)+Ql(N2, 2)*X(N2) 

+Q 1 (N2 , 3 ) *X(N7 ) +Q 1 (N2 , N3 ) *EXP ( EK-EK/ 

X  ( N  2  )  )*X(2)/(l.  +DKV  1*X(  2  )  ) 

220 

DO  890  I=N7 ,N6 

221 

L2=I-N2 

222 

L3=I-N3 

223 

I 

i— i 

II 

224 

L5=I-N5 

225 

890 

G ( I ) =Q 1 ( I , L2 ) *X( 1-2 ) +Q 1 ( I , L3 ) *X( 1-1 ) + 
Q1(I,L4)*X(I)+Q1(I,L5)*X(I+1)+Q1(I,N3) 
*EXP(EK-EK/X(l)  )  *X(  L4  )/(  1 .  +DKV1*X( L4  )  ) 

226 

G(Nl)=Ql(Nl,N4-2)*X(Nl-2)+Ql(Nl,N4-l) 
*X(N1-1)+Q1(N1,N4)*X(N1)+Q1(N1,N2)+ 
Ql(Nl ,N3)*EXP(EK-EK/x(N1) )*X(N4)/ 

( 1 . +DKV 1*X(N4 ) ) 

227 

C 

SUMMING  THE  NEW  RESIDUALS 

228 

G1=0 . 0 

229 

DO  599  1=1, N1 

230 

G2=ABS ( G ( I ) ) 

231 

599 

G1=G2+G1 

232 

C4(I1)=G1 

233 

C 

CHECKING  FOR  CONVERGENCE 

234 

C 

COMPARING  THE  RESIDUALS  SUM 

235 

IF(C4(I1).GE.C4(I1-1))G0  TO  68 

' 


■ 

■ 


237 

238 

239 

240 

241 

242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

265 

266 
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GO  TO  67 

C  CALCULATE  NEW  INCREMENT  IF  NECESSARY 

68  12=12+1 

13=13+1 

IF(I3.LT.I2)G0  TO  888 
Dl(l2)=Dl(l2-l)/2 
DO  171  1=1, N1 
171  X(I)=F(I) 

GC  A  iOGO  TO  100 
888  15=15+1 

'.-15  Bl  =  I5 

D1(I2)=B+C6/B1 
DO  175  1=1, Ni 
175  X( I ) =F ( I ) 

GO  TO  100 
67  v=.0005 

C  CHECKING  THE  SIZE  OF  EACH  RESIDUAL 

DO  555  1=1, Ni 
IF(ABS(G(I) ) .GT.V)GOTO  40 
555  CONTINUE 

B=D1 ( 12 ) 

15=1 

IF(Dl(l2) . GE . W ) GO  TO  150 
CONTINUE 
DO  173  1=1, Ni 
173  F ( I ) =X( I ) 

12=12+1 

D1  ( 12 ) =D1 ( 12-1 ) *2 . 

C6=ABS (Dl( 12-2 )-Dl( 12-1) ) 

GO  TO  100 
150  CONTINUE 


' 


. 


■ 


. 


' 


NUMB 

267 

268 

269 

270 

271 

272 

273 

274 

275 

276 

277 

278 

279 

280 

281 

282 

283 

284 

285 

286 

287 

288 

289 

290 
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DO  176  1=1, N1 
176  F(T)=X(I) 

C  EVALUATION  OF  LOCAL  EFFECTIVENESS 

FACTORS 

RO=DKV*PX/(  1 .  +DICV  1*PX) 

DO  73  1=1, N4 
NN2=N3-1+I 

73  R(I)=DKV*EXP(EK-EK/X(NN2) )*X(I)/  ' 

( 1  .  +DICV  1*X(  I )  ) 

DO  74  1=1, N4 

74  EFF(I)=R(I)/R0 

C  EVALUATION  OF  OVERALL  EFFECTIVENESS 

FACTOR 

dimensiona3(5,6) ,  A2 ( 5 , 6 ) , ai ( 5 , 6 ) , 

Yi(5),Y2(5),Y3(5) 

DO  800  J=l,6 

800  A3(l, J)=1.000 
DO  801  1=2,5 

801  A3(I,1)=0.0000 
DO  802  1=2,5 

802  A3 ( I , 2 ) =1 . /( 4 . ** ( 1-1 ) ) 

DO  803  1=2,5 

803  A3(I,3)=1./(2.**(I-1) ) 

DO  804  1=2,5 

804  A3(I,4)=3.**(I-1)/(4.**(I-1) ) 

DO  805  1=2,5 

805  A3(I,5)=1.000 
DO  806  1=2,5 


. 


NUMB 

291 

292 

293 

294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 

305 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315 

316 

317 

318 
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99=1 

806 

A3(I,6)=1./B9 

DO  807  1=1,5 

DO  807  J=1 , 6 

807 

A2(I, J)=A3(I, J) 

DO  808  1=1,5 

DO  808  J=1 , 6 

808 

A 1 ( I , J ) =A3 ( I , J ) 

DO  809  1=2,5 

809 

A2(I,2)=1./(7.**(I-1) ) 

DO  810  1=2,5 

810 

A2(l,3)=3.**(l-i)/(7.**(l-i)) 
DO  811  1=2,5 

811 

A2(I,4)=5.**(I-1)/(7.**(I-1)) 

DO  812  1=2,5 

812 

Al(l,2)=2.**(l-l)/(7 .**(1-1) ) 
DO  813  1=2,5 

813 

A1(I,3)=4.**(I-1)/(7.**(I-1)) 
DO  814  1=2,5 

814 

A1(I,4)=6.**(I-1)/(7.**(I-1) ) 

C 

SOLVE  THE  MATRIX  A3,A2,A1 

DO  815  1=1,5 

DO  815  J=1 , 6 

815 

A ( I , J ) =A3 ( I , J ) 

N=5 

110  =  2 

19  =  2 

. 

GO  TO  19 

' 


. 


•  r  ct-  .  ■  .  i 

. 
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319 

816 

CONTINUE 

320 

DO  817  1=1,5 

321 

817 

Y3(I)=X(I)*7./DN1 

322 

DO  818  1=1,5 

323 

DO  818  J=1 , 6 

324 

818 

A ( I j  J ) =A2 ( I , J ) 

325 

N=5 

326 

110=3 

327 

w 

00 

II 

328 

GO  TO  19 

329 

819 

CONTINUE 

330 

DO  820  1=1,5 

331 

820 

Y2 ( I ) =X( I ) *8 ,/DNl 

332 

DO  821  1=1,5 

333 

DO  821  J=1 , 6 

334 

821 

A ( I , J ) =A1 ( I , J ) 

335 

N=5 

336 

110=4 

337 

17  =  4 

338 

GO  TO  19 

339 

822 

CONTINUE 

340 

• 

DO  823  1=1,5 

341 

823 

yi(i)=x(i)*7./dni 

342 

DN 10=0.0 

343 

1  =  2 

344 

J  =  1 

345 

t 

ii 

M- 

346 

GO  TO  831 

' 

, 


^cv. 


• 

• 
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347 

824 

1  =  1  +  1 

348 

♦ 

‘J=J  +  1 

349 

K=K+2 

350 

831 

DN9=K 

351 

DN1 1=Y1 ( I ) *EFF ( J ) *DN9**2 

352 

DN10=DN10+DN 1 1 

353 

IF ( I . EQ . 4 ) GO  TO  830 

354 

GO  TO  824 

355 

830 

DN12=0 . 00 

356 

B 1 3  = Y l(5) + Y  2(l) 

357 

NN4=N4-7 

358 

NN5=N1-15 

359 

M 

1! 

360 

J=7 

361 

GO  TO  832 

362 

825 

1=1  +  4 

363 

J= J+8 

364 

832 

DN9=J 

365 

DN13=B13*EFF(l)*DN9**2 

366 

DN 1 2  =DN 1 3 +DN 1 2 

367 

DN25=Nl-7 

368 

, 

IF ( I . EQ .NN4 ) GO  TO  833 

369 

833 

DN14= ( Y2 ( 5 ) +Y3 ( 1 ) ) *EFF (NN4-3 ) *DN25**2 

370 

DNl6=0.000 

371 

NN6=N4-6 

372 

NN7=N4-5 

373 

• NN8=N4-4 

374 

NN9=N1-13 

375 

NN10=N1-12 

. 


. 


NUMB 

376 

377 

378 

379 

380 

381 

382 

383 

384 

385 

386 

387 

388 

389 

390 

391 

392 

393 

394 

395 

396 

397 

398 

399 

400 

401 

402 

403 
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NNli=Nl-ll 

J=5 

K=9 

GO  TO  834 

826  J= J+4 
K=K+8  ' 

834  DN9=K 
DN15=Y2(2)*EFF(J)*DN9**2 
DNl6=DNl5+DNl6 

IF( J.EQ .NN6)G0  TO  835 
GO  TO  826 

835  DN15=0.00 
J  =  6 

K=10 

GO  TO  836 

827  J=J+4 
K=K+8 

836  DN9=K 
DN17=Y2(3)*EFF(J)*DN9**2 
IF(J.EQ .NN7)G0  TO  837 

GO  TO  827 

837  DN13=0 . 000 

J=7 

K=ll 

GO  TO  838 

828  J=J+4 
K=K+4 

838  DN9=K 
DN18=Y2(4)*EFF(l )*DN9**2 


. 


UMB 

405 

406 

407 

408 

409 

410 

411 

412 

413 

414 

415 

4l6 

417 

418 

419 

420 

421 

422 

423 

424 

425 

426 

427 

428 

429 

430 

431 


(87) 


STATEMENT  STATEMENT 

'  NUMBER 


DN 1 3  =DN 1 3  +DN 1 8 

IF(J.EQ .NNll)G0  TO  839 

GO  TO  828 

DN20=Nl-5 

DN2i=Nl-3 

DN22=N1- 1 

DN19=Y3 ( 2 ) *EFF (N9 ) *DN20**2+Y3 ( 3 ) *EFF ( 

% 

N5)*DN21**2+Y3(4)*EFF(N4)*DN22**2+Y3(5) 

EFFET=DN10+DN12+DNl6+DN17+DN13+DNl9 

DN26=N1 

EFFO= ( 3 . /DN26**2 ) *EFFET 

THMOD=EFFO*DKV*RG*TS/(DEX*  ( 1 .  +DKV1*PX)  ) 
BETA=CXS*DELH*DEX/(DALP*TS ) 

GAMMA=EK 

CF=PX*DKV1 

93  FORMAT ( IX, 2 1HGAMMA, BETA, TH. MOD. ,CF) 

WRITE  (6,93) 

2  FORMAT(lX,9F14.5) 

WRITE (6,2) GAMMA , BETA , THMOD , CF 
91  F0RM4T ( IX, 27HL0CAL  EFFECTIVENESS  FACTORS) 

WRITE (6, 2) (EFF(I) , I=1,N4) 

105  FORMAT ( IX, 37HEFFECTIVENESS  FACTOR  AND, 

DEX , DALP , DKV ) 

WRITE (6,2) EFFO , DEX, DALP , DKV 

6  FORMAT ( IX, 9HRESIDUALS) 

WRITE (6, 6) 

WRITE (6, 2) (G(l) , I=1,N1) 

7  FORMAT ( IX, 29HC0NCENTRAT ION  AND  TEMPERATURE) 
WRITE (6, 7 ) 

WRITE (6, 2) (X(l) , I = 1 , N 1 ) 


' 


. 


■ 


(88) 


CARD 

NUMBER 

STATEMENT 

NUMBER 

STATEMENT 

433 

126 

FORMAT ( IX, 50HSURFACE  TEMP, PARTIAL  PRESS 
, HEAT  OF  REAX, DELTA  R  SQ ) 

434 

WRITE (6, 126) 

435 

WRITE ( 6 , 2 ) TS , PX , DELH , DELR 

436 

STOP 

437 

END 

***  The  following  corrections  are  required  to  make  the 
program  work  (Program  l). 

Between  card  numbers  207  and  208 
IF(I10.EQ.I9)G0  TO  816 
IF(li0.EQ.I8)G0  TO  825 
Between  card  numbers  368  and  369 


GO  TO  825 


■ 


■ 


Program  2 


(«9) 


CARD  STATEMENT  STATEMENT 

NUMBER'  '  NUMBER 


24 

P(l, l)=Z*(-i./DELR) 

25 

P ( 1 »  2 )  =-P ( 1 , 1 ) 

26 

P(2 , 1 ) =Z*DELR 

27 

P(2,2)=Z*(-2./DELR) 

28 

P(2 , 3)=Z*(  1  ./DELR) 

32 

NO  EQUIVALENT  CARD... 

33 

P(l, i-i)=z/delr 

34 

P ( I , I )=Z*(-2./DELR) 

35 

p(i , i+i)=z/delr 

37 

P(N4,N9)=Z*(-1./(5.*DELR) ) 

38 

P(N4,N5)=Z*(lO./(5 .*delr) ) 

39 

P(N4,N4)=Z*(-25./(5.*DELR) ) 

44 

P(N4,N2)=Z*(l6./(5.*DELR) )*px 

45 

P(Nl,N2)=TS*DALP*(l6./(5.*DELR) ) 

**  Cards  277 

to  410  inclusive  have  no  equivalent  in  Program  2 

411 

EFFET=0 . 000 

412 

D0  75  1=1, N4 

413 

75  EFFET=EFF ( I ) +EFFET  1 

414 

EFF0=2/DN1*EFFET 

. 


( * ..  I 


■ 


. 

. 


; 


. 


